text:
# Deterministic computations whose history is independent of the order of asynchronous updating

## 1. Introduction

Consider a set $``$ of processors (sites) in which each site $`x`$ has a set $`\mathrm{SS}`$ of possible states (also called “local states”). An arbitrary function $`\xi \mathrm{SS}^{}`$ is called a space-configuration, or simply “configuration”, or “global state”. The value $`\xi (x)`$ is the state of site $`x`$ in $`\xi `$. A function $`N:2^{}`$ will be called a neighborhood function assigning to each $`x`$, a set $`N(x)`$ called the neighborhood of $`x`$. A function $`f:\mathrm{SS}^{}\mathrm{SS}^{}`$ is called a transition function if $`f(\xi )(x)`$ depends only on $`\xi N(x)`$, i.e.

$$\xi _1N(x)=\xi _2N(x)f(\xi _1)(x)=f(\xi _2)(x).$$

The transition function determines a possible “next” configuration from the “current” one. The 4-tuple

(1.1)
$$𝔸=(,\mathrm{SS},N,f)$$

will be called an automaton (not necessarily a finite one). If all sets $`N(x)`$ are finite then the system is called local. Note that locality is actually a property of $`f`$ itself: it says that for each $`x`$ a finite $`N(x)`$ can be chosen such that $`f(\xi )(x)`$ depends only on $`\xi N(x)`$. Let $`_+=[0,\mathrm{})`$.

###### Example 1.1 (Cellular automata).

1. On the set of integers: Let $`=`$, $`N(x)=\{x1,x,x+1\}`$. Suppose that there is a transition function $`g(x,y,z)`$ and each site $`x`$ has a value $`\xi (x)\mathrm{SS}`$. Now the result of transition at site $`x`$ is

$$f(\xi )(x)=g(\xi (x1),\xi (x),\xi (x+1)).$$

In this example, the transition function depends only on the sequence of values of $`\xi N(x)`$, i.e. it is homogenous. The present paper will not exploit any consequences of homogeneity.
2. On the set of natural numbers, with “free boundary condition”: Let $`=_+`$, $`N(x)=\{x1,x,x+1\}`$ for $`x>0`$ and $`\{0,1\}`$ for $`x=0`$. Suppose that there are transition functions $`g(x,y,z)`$, $`g_0(x,y)`$. Now the result of transition at site $`x`$ is $`g(\xi (x1),\xi (x),\xi (x+1))`$ for $`x>0`$ and $`g_0(\xi (0),\xi (1))`$ for $`x=0`$.

Let us fix an automaton $`𝔸`$ as in (1.1). An arbitrary function $`\eta :\times _+\mathrm{SS}`$ is called a space-time configuration. Such a space-time configuration can also be viewed as a sequence $`\eta :_+\mathrm{SS}^{}`$ of space-configurations. We will say that a space-time configuration $`\eta `$ is a synchronous trajectory if for all $`x,t`$ we have $`\eta (,t+1)=f(\eta (,t))`$. In other words,

(1.2)
$$\eta (x,t+1)=f(\eta (,t))(x),$$

i.e. in $`\eta `$, each site is “updated” every time by the function $`f`$ (though the update may not change the state). We are interested in situations when at any one time, only the values of some of the sites are updated. We will say $`\eta `$ is an asynchronous trajectory if (1.2) holds for all $`x,t`$ such that $`\eta (x,t+1)\eta (x,t)`$: i.e. if each site in $`\eta `$ at each time is either updated or left unchanged. From now on, when we speak of a “trajectory” without qualification, this will mean an asynchronous trajectory. Let the update set

$$U(t,\eta )$$

be the set of sites $`x`$ with $`\eta (x,t+1)\eta (x,t)`$. The initial configuration and the update sets $`U(t,\eta )`$ determine $`\eta `$. For any set $`A`$, let

$$\chi (x,A)=\{\begin{array}{cc}1\hfill & \text{if }xA\text{,}\hfill \\ 0\hfill & \text{otherwise.}\hfill \end{array}$$

For a space-time configuration $`\eta (x,t)`$ we define the function $`\tau (x,t)=\tau (x,t,\eta )`$ as follows:

$`\tau (x,0)`$ $`=0,`$
$`\tau (x,t+1)`$ $`=\tau (x,t)+\chi (x,U(t,\eta )).`$
We can call $`\tau (x,t)`$ the effective age of site $`x`$ in the space-time configuration $`\eta `$ at time $`t`$: this is the number of effective updatings that $`x`$ underwent until time $`t`$. Given an initial configuration $`\xi `$, we say that $`f`$ (and thus $`𝔸`$) has invariant histories on $`\xi `$ if there is a function $`\zeta (x,u)=\zeta (x,u,\xi )`$ such that for all asynchronous trajectories $`\eta (x,t)`$ with $`\eta (,0)=\xi `$ we have

(1.3)
$$\eta (x,t)=\zeta (x,\tau (x,t,\eta ),\xi ).$$

This means that after eliminating repetitions, the sequence $`\zeta (x,0),\zeta (x,1)`$, $`\mathrm{}`$ of values that a site $`x`$ will go through during some space-time configuration, does not depend on the update sets, only on the initial configuration (except that the sequence may be finite if there is not an infinite number of successful updates). The update sets influence only the delays in going through this sequence. We say that an automaton has invariant histories if it has such on all initial configurations.

###### Remark 1.2.

The sequence $`\zeta (x,0),\zeta (x,1),\mathrm{}`$ is a sequence of local states but $`\zeta (,n)`$ is not a space-configuration (global state) that appears at any time in a typical asynchronous trajectory. ∎

###### Theorem 1.3.

If $`𝔸`$ is a one-dimensional cellular automaton with state space $`\mathrm{SS}=\{0,\mathrm{},n1\}`$ for some natural number $`n`$, then it is undecidable whether it has invariant histories.

The theorem shows that some more condition is needed if we want the invariant history property to become decidable. For us, this condition will be monotonicity. The set of free sites $`x`$ in a configuration $`\xi `$ is defined by

$$L(\xi )=\{x:f(\xi )(x)\xi (x)\}.$$

For a space-time configuration $`\eta `$, let

$$L(t,\eta )=L(\eta (,t)).$$

For a configuration $`\xi `$ and a set $`E`$ of sites, let

$`f(\xi ,E)(x)`$ $`=\{\begin{array}{cc}f(\xi )(x)\hfill & \text{if }xE\hfill \\ \xi (x)\hfill & \text{otherwise.}\hfill \end{array}`$
$`f(\xi ,E,F)`$ $`=f(f(\xi ,E),F).`$
With this notation, we have $`f(\xi )=f(\xi ,)=f(\xi ,L(\xi ))`$. Now we can express the condition that $`\eta `$ is an asynchronous trajectory by saying that for every $`t`$ there is a set $`U`$ with

(1.4)
$$\eta (,t+1)=f(\eta (,t),U),$$

and the condition that $`\eta `$ is synchronous by requiring $`U(t,\eta )=L(t,\eta )`$ for each $`t`$. We call a transition rule $`f`$ monotonic if $`L(t,\eta )U(t,\eta )L(t+1,\eta )`$, i.e. updating a site cannot take away the freedom of other sites. We call a transition rule $`f`$ (and thus the automaton $`𝔸`$) commutative if for all configurations $`\xi `$ and all disjoint sets of sites $`A,BL(\xi )`$ we have

(1.5)
$$f(\xi ,A,B)=f(\xi ,AB).$$

We call $`f`$ locally commutative when this property is required just for the special case where $`A,B`$ are one-element sets. The following fact is easy to see but we give the proof for completeness.

###### Lemma 1.4.

If $`f`$ is local then its local commutativity implies commutativity.

###### Proof.

Let us first show

(1.6)
$$f(\xi ,\{x_1\},\mathrm{},\{x_n\})=f(\xi ,\{x_1,\mathrm{},x_n\}).$$

Local commutativity implies for each $`k`$,

$$\xi ^{}=f(\xi ,\{x_1\},\mathrm{},\{x_n\})=f(\xi ,\{x_k\},\{x_1\},\mathrm{},\{x_{k1}\},\{x_{k+1}\},\mathrm{},\{x_n\}).$$

Therefore $`\xi ^{}(x_k)=f(\xi ,\{x_1,\mathrm{},x_n\})(x_k)`$. Now, let us show

(1.7)
$$f(\xi ,\{x_1,\mathrm{},x_n\},\{y\})=f(\xi ,\{x_1,\mathrm{},x_n,y\}).$$

Using (1.6), we have $`f(\xi ,\{x_1\},\mathrm{},\{x_n\})=f(\xi ,\{x_1,\mathrm{},x_n\})`$, hence $`f(\xi ,\{x_1,\mathrm{},x_n\},\{y\})=f(\xi ,\{x_1\},\mathrm{},\{x_n\},\{y\})`$. Using (1.6) again concludes the proof.

Let us return to the general case. Obviously, it is sufficient to check (1.5) for sites $`yB`$. Clearly, $`f(\xi ,A,B)(y)=f(\xi ,N(y)A,\{y\})`$. The latter is $`f(\xi ,(N(y)A)\{y\})`$ according to (1.7). ∎

###### Remarks 1.5.

1. For the cellular automaton example above, local commutativity is equivalent to saying that if $`g(r_0,r_1,r_2)r_1`$ and $`g(r_1,r_2,r_3)r_2`$ then

$`g(g(r_0,r_1,r_2),r_2,r_3)`$ $`=g(r_1,r_2,r_3)`$
$`g(r_0,r_1,g(r_1,r_2,r_3))`$ $`=g(r_0,r_1,r_2).`$
2. If $`f`$ is not local then local commutativity does not always imply commutativity. For an example, let $`=\{0,1\}`$, $`=`$, $`N(x)=`$, and let

$$f(\xi )(x)=\{\begin{array}{cc}1\hfill & \text{if }\xi (y)=0\text{ for all but finitely many }y\text{,}\hfill \\ 0\hfill & \text{otherwise.}\hfill \end{array}$$

Now $`f`$ is obviously locally commutative. On the other hand, let $`\xi _0(x)=0`$ for all $`x`$, and let Then $`f(\xi _0,)(1)=1`$ and $`f(\xi _0,_+,_+)(1)=0`$.

###### Theorem 1.6.

A transition function is commutative if and only if it is monotonic and has invariant histories.

In Theorem 3.1 below, we will give a known simple example of a commutative transition function. For that example, the theorem can be proved much easier.

Theorem 1.6 can be derived from results e.g. in . However, I do not find it worth introducing all the concepts needed for the derivation: the simplicity of the condition in the present context probably justifies a self-contained proof.

## 2. Commutativity implies invariant histories

###### Lemma 2.1.

Suppose that $`f`$ has invariant histories and is monotonic: then it is commutative.

###### Proof.

Let $`U_1(0)=U_2(1)=\{x\}`$, $`U_1(1)=U_2(0)=\{y\}`$, and $`U_1(t+2)=U_2(t+2)`$. This defines $`\eta _1`$ and $`\eta _2`$ from initial configuration $`\xi `$ by $`U_1,U_2`$ as usual. By monotonicity, $`\eta _1(y,1)\eta _1(y,2)`$ and $`\eta _2(x,1)\eta _2(x,2)`$, so $`\tau `$’s values satisfy

$$\tau (x,2,\eta _1)=\underset{t=0}{\overset{1}{}}\chi (w,U_1(t,\eta _1))$$

which is 1 if $`w\{x,y\}`$ and 0 otherwise. The same value is obtained for $`\tau (x,2,\eta _2)`$. By invariant histories, there is a $`\zeta `$ such that

$$\eta _1(w,2)=\zeta (w,\tau (w,2,\eta _1))=\zeta (w,\tau (w,2,\eta _2))=\eta _2(w,2)$$

and

$$f(\xi ,x,y)=f(\xi ,U_1(0),U_1(1))=\eta _1(w,2)=\eta _2(w,2)=f(\xi ,U_2(0),U_2(1))=f(\xi ,\{y\},\{x\}).$$

Thus, $`f`$ is commutative. ∎

What remains to prove after Lemma 2.1 is that commutativity implies monotonicity and invariant histories.

###### Lemma 2.2.

If $`f`$ is commutative then it is monotonic.

###### Proof.

By Lemma 1.4, $`f(\xi ,U(t,\eta ),L(t,\eta )U(t,\eta ))=f(\xi ,L(r,\eta ))`$. Therefore $`L(t,\eta )U(t,\eta )L(t,\eta )`$ implies that $`f`$ is monotonic. ∎

We say for two asynchronous trajectories $`\eta _0,\eta _1`$ with the same initial configuration that $`\eta _1`$ dominates $`\eta _0`$ until time $`u`$ if the following conditions hold:

1. $`\tau (,t,\eta _0)\tau (,t,\eta _1)`$ for all $`tu`$
2. for all $`t_0,t_1u`$, if $`\tau (x,t_0,\eta _0)=\tau (x,t_1,\eta _1)`$ then $`\eta _0(x,t_0)=\eta _1(x,t_1)`$.

When $`\eta _1`$ dominates $`\eta _0`$ up to time $`u`$ for all $`u`$ then we simply say that $`\eta _1`$ dominates $`\eta _0`$. This domination is, of course, a transitive relation. If the rule has invariant histories then condition (a) implies (b), but otherwise this may not be the case.

###### Proof of Theorem 1.6.

Let $`f`$ be a commutative transition rule. It remains to prove that it has invariant histories.

Let $`\eta `$ be an asynchronous trajectory and $`A_0L(0,\eta )U(0,\eta )`$. Then there is an asynchronous trajectory $`\eta ^{}`$ dominating $`\eta `$ with initial configuration $`\eta (,0)`$, such that $`U(0,\eta ^{})=U(0,\eta )A_0`$.

Proof. Let $`\xi _0=\eta (,0)`$. We show how to build, for each $`u`$, a trajectory $`\eta ^{}`$ with the given properties that dominates $`\eta `$ up to time $`u`$. When $`u\mathrm{}`$ then $`\eta ^{}`$ will converge to a trajectory with the same properties that dominates $`\eta `$. For $`u=0`$ we can choose $`\eta ^{}(,0)=\eta (,0)`$. We assume that $`\eta ^{}`$ can be constructed for all $`v<u`$ and prove it for $`u`$. Let $`\xi _1=\eta (,1)`$, and $`A_1=A_0U(1,\eta )`$. Let the trajectory $`\eta _1`$ be defined by $`\eta _1(x,t)=\eta (x,t+1)`$. The inductive assumption gives a trajectory $`\eta _1^{}`$ with initial configuration $`\xi _1`$ dominating $`\eta _1`$, with

(2.1)
$$U(0,\eta _1^{})=A_1U(0,\eta _1).$$

Using this trajectory, we define, for $`t>0`$:

$$\eta ^{}(,t)=\{\begin{array}{cc}f(\xi _0,A_0U(0,\eta ))\hfill & \text{if }t=1\text{,}\hfill \\ \eta _1^{}(,t1)\hfill & \text{otherwise.}\hfill \end{array}$$

$`\eta ^{}`$ is an asynchronous trajectory.

Proof. Let us show that $`\eta ^{}`$ satisfies (1.4). This holds by definition for $`t=0`$ and $`t>1`$. Let us show that it also holds for $`t=1`$ with $`U=U(1,\eta )A_0`$. We have

(2.2) $`\eta ^{}(,2)`$ $`=\eta _1^{}(,1)`$ $`\text{by def.},`$
$`=f(\xi _1,A_1U(0,\eta _1))`$ $`\text{by (}\text{2.1}\text{)},`$
$`=f(\xi _1,A_1U(1,\eta ))`$ $`\text{by def. of }\eta _1,`$
$`=f(\xi _1,(A_0U(1,\eta ))U(1,\eta ))`$ $`\text{by def. of }A_1,`$
$`=f(\xi _1,A_0(U(1,\eta )A_0))`$
$`=f(\xi _0,U(0,\eta ),A_0,U(1,\eta )A_0)`$ $`\text{by def. of }\xi _1\text{ and commutativity},`$
$`=f(\eta ^{}(,1),U(1,\eta )A_0).`$
For domination, we must check two properties.

We have $`\tau (x,t,\eta )\tau (x,t,\eta ^{})`$.

Proof. By the definition of $`\tau `$, for $`t>0`$,

$$\tau (x,t,\eta )=\{\begin{array}{cc}\chi (x,U(0,\eta ))\hfill & \text{if }t=1\text{,}\hfill \\ \tau (x,1,\eta )+\tau (x,t1,\eta _1)\hfill & \text{if }t>1\text{.}\hfill \end{array}$$

By the definition of $`\eta _1^{},\eta `$, for $`t>0`$, using (2.2), we have

(2.3)
$$\tau (x,1,\eta _1^{})=\chi (x,A_1U(1,\eta ))=\chi (x,A_0U(1,\eta )).$$

Further,

(2.4)
$$\tau (x,t,\eta ^{})=\{\begin{array}{cc}\chi (x,A_0U(0,\eta ))\hfill & \text{if }t=1\text{,}\hfill \\ \tau (x,1,\eta ^{})+\chi (x,U(1,\eta )A_0)\hfill & \text{if }t=2\text{,}\hfill \\ \tau (x,2,\eta ^{})+\tau (x,t1,\eta _1^{})\tau (x,1,\eta _1^{})\hfill & \text{if }t>2\text{.}\hfill \end{array}$$

By the above definition,

$`\tau (x,1,\eta ^{})`$ $`=\tau (x,1,\eta )+\chi (x,A_0),`$
$`\tau (x,2,\eta ^{})`$ $`=\tau (x,1,\eta )+\chi (x,A_0)+\chi (x,U(1,\eta )A_0)`$
$`=\tau (x,1,\eta )+\chi (x,A_0U(1,\eta ))\tau (x,1,\eta )+\chi (x,U(1,\eta ))`$
$`=\tau (x,2,\eta ).`$
Also, from here and (2.3),

(2.5)
$$\tau (x,2,\eta ^{})=\tau (x,1,\eta )+\chi (x,A_0U(1,\eta ))=\tau (x,1,\eta )+\tau (x,1,\eta _1^{}).$$

By domination, $`\tau (x,t1,\eta _1^{})\tau (x,t1,\eta _1)`$ and hence for all $`t2`$, we have, combining (2.4) with (2.5),

(2.6) $`\tau (x,t,\eta ^{})`$ $`=\tau (x,1,\eta )+\tau (x,t1,\eta _1^{})`$
$`\tau (x,1,\eta )+\tau (x,t1,\eta _1)=\tau (x,t,\eta ).`$
If $`\tau (x,s,\eta )=\tau (x,s^{},\eta ^{})`$ then $`\eta (x,s)=\eta ^{}(x,s^{})`$.

Proof. If $`\tau (x,s,\eta )=0`$ then clearly $`\eta ^{}(x,s)=\eta ^{}(x,s^{})`$ since this means that in both processes, no progress has been made in $`x`$ from the initial configuration. Assume therefore that $`\tau (x,s,\eta )>0`$ and hence $`s,s^{}>0`$.

Assume $`s^{}=1`$. Then $`\tau (x,s,\eta )=\tau (x,1,\eta ^{})=1`$ and hence $`xA_0U(0,\eta )`$. If $`xU(0,\eta )`$ then $`s=1`$ and hence the same transition that gives $`\eta ^{}(x,1)`$ also gives $`\eta (x,1)`$. Otherwise $`s>1`$ hence $`\tau (x,s1,\eta _1)=1`$. Also, $`xA_0U(0,\eta _1^{})`$, hence $`\tau (x,1,\eta _1^{})=1`$. The inductive assumption implies $`\eta _1^{}(x,1)=\eta _1(x,s1)=\eta (x,s)`$. On the other hand, (2.2) and $`xU(0,\eta )`$ implies $`\eta _1^{}(x,1)=\eta ^{}(x,1)`$ which concludes this case.

Assume now $`s^{}>1`$. Since $`\eta (x,t)`$ changes if and only if $`\tau (x,t)`$ does we can assume that $`xU(s,\eta )`$ since otherwise we can decrease $`s`$ without changing $`\eta (x,s)`$. The same is true for $`s^{}`$. Under these assumptions we have $`ss^{}`$. By (2.6),

$$\tau (x,s^{},\eta ^{})=\tau (x,1,\eta )+\tau (x,s^{}1,\eta _1^{}).$$

We assumed this to be equal to $`\tau (x,s,\eta )=\tau (x,1,\eta )+\tau (x,s1,\eta _1)`$. Hence $`\tau (x,s^{}1,\eta _1^{})=\tau (x,s1,\eta _1)`$. Also $`\eta (x,s)=\eta _1(x,s1)`$, $`\eta ^{}(x,s^{})=\eta _1^{}(x,s^{}1)`$, and hence the inductive assumption implies the statement.

Let $`\eta `$ be a trajectory. Then the synchronous trajectory with initial configuration $`\eta (,0)`$ dominates $`\eta `$.

Proof. Let $`A_0=L(0,\eta )U(0,\eta )`$. By 2 above, there is a trajectory $`\eta ^{}`$ with initial configuration $`\eta (,0)`$ dominating $`\eta `$ such that $`U(0,\eta ^{})=U(0,\eta )A_0=L(0,\eta )`$. This just means that $`\eta ^{}`$ is a synchronous trajectory up to time 1. Continuing the application of 2, we can dominate $`\eta `$ by a synchronous trajectory $`\eta ^{\prime \prime }`$ up to time 2, etc.

Now we can conclude the proof of the theorem as follows. Let $`\eta `$ be a trajectory with initial configuration $`\xi `$ and let $`\eta ^{}`$ be the synchronous trajectory with the same initial configuration. Let us define

$`\sigma (x,s,\xi )`$ $`=\mathrm{min}\{t:\tau (x,t,\eta ^{})=s\},`$
$`\zeta (x,s,\xi )`$ $`=\eta ^{}(x,\sigma (x,s)).`$
To prove (1.3), note that due to domination, $`\tau (x,t,\eta )\tau (x,t,\eta ^{})`$ and hence for every $`x,y,t`$ there is a $`t^{}t`$ with $`\tau (x,t,\eta )=\tau (x,t^{}\eta ^{})`$. Let $`t^{}`$ be the first such: $`t^{}=\sigma (s,\tau (x,t,\eta ))`$. By domination, $`\eta (x,t)=\eta ^{}(x,t^{})=\zeta (x,t)`$. ∎

## 3. A rich example of commutative transitions

Let us show the known result that every transition function can be embedded into a commutative one. We will use the following notation:

$$bamodm$$

is the integer $`x`$ with $`xb(modm)`$ and $`m/2<xm/2`$.

###### Theorem 3.1.

Let $`𝔸_1=(,\mathrm{SS}_1,N,f_1)`$ be an arbitrary local (not necessarily commutative) automaton $`N(x)`$. Then there is an automaton $`𝔸_2=(,\mathrm{SS}_1\times R,N,f_2)`$, where for $`s\mathrm{SS}_1\times R`$ we write $`s=(s.\text{F},s.\text{G})`$, with the following property.

Let $`\xi _1`$ be an arbitrary configuration of $`f_1`$ and let $`\xi _2`$ be a configuration of $`f_2`$ such that for all $`x`$ we have $`\xi _2(x)=(\xi _1(x),0\mathrm{}0)`$. Then for the synchronous trajectory $`\eta _2`$ of $`f_2`$, with initial configuration $`\xi _2`$, the space-time configuration $`(x,t)\eta _2(x,t).\text{F}`$ is a synchronous trajectory of $`f_1`$. Moreover, in this trajectory, the state of each cell changes in each step.

In other words, as long as we update synchronously the rule $`f_2`$ behaves in its field F just like the arbitrary rule $`f_1`$. But $`f_2`$ has invariant histories, so it is much more robust.

###### Proof.

Let $`\mathrm{SS}_2=\mathrm{SS}_1^2\times \{0,1,2\}`$. The three components of each state $`s`$ of $`\mathrm{SS}_2`$ will be written as

$$s.\text{Cur},s.\text{Prev}\mathrm{SS}_1,s.\text{Age}\{0,1,2\}.$$

The statement of the theorem will obtain by $`s.\text{F}=s.\text{Cur}`$, $`s.\text{G}=(s.\text{Prev},s.\text{Age})`$. The field $`\text{Age}\{0,1,2\}`$ will be used to keep track of the time of the simulated cells mod 3, while Prev holds the value of Cur for the previous value of Age.

Let us define $`s^{}=f_2(\xi )(x)`$. If there is a $`yN(x)`$ such that $`(\xi (y).\text{Age}\xi (x).\text{Age})amod3<0`$ (i.e. some neighbor lags behind) then $`s^{}=\xi (x)`$ i.e. there is no effect. Otherwise, let $`\sigma (y)`$ be $`\xi (y).\text{Cur}`$ if $`\xi (y).\text{Age}=\xi (x).\text{Age}`$, and $`\xi (y).\text{Prev}`$ otherwise.

$`s^{}.\text{Cur}`$ $`=f_1(\sigma )(x),`$
$`s^{}.\text{Prev}`$ $`=\xi (x).\text{Cur},`$
$`s^{}.\text{Age}`$ $`=\xi (x).\text{Age}+1mod3.`$
Thus, we use the Cur and Prev fields of the neighbors according to their meaning and update the three fields according to their meaning. It is easy to check that this transition rule simulates $`f_1`$ in the Cur field if we start it by putting 0 into all other fields.

Let us check that $`f_2`$ is locally commutative. If two neighbors $`x,y`$ are both are allowed to update then neither of them is behind the other modulo 3, hence they both have the same Age field. Suppose that $`x`$ updates before $`y`$. In this case, $`x`$ will use the the Cur field of $`y`$ for updating and put its own Cur field into Prev. Next, since now $`x`$ is “ahead” according to Age, cell $`y`$ will use the Prev field of $`x`$ for updating: this was the Cur field of before. Therefore the effect of consecutive updating is the same as that of simultaneous updating. ∎

The commutative medium of the above proof is also called the “marching soldiers” scheme since its handling of the Age field reminds one of a chain of soldiers marching ahead in which two neighbors do not want to be separated by more than one step. It is shown in  that if the update times obey a Poisson process then the average computation time of this simulation within a constant factor of the computation time of the synchronous computation.

###### Remark 3.2.

In typical cases of asynchronous computation, there are more efficient ways to build a commutative rule than to store the whole previous state in the Prev field. Indeed, the transition function typically does not use the complete state of cells in $`N(x)`$. Rather, the cells only “communicate” in the sense that there is a message field and the next state of $`x`$ depends only on this field of the neighbor cells. In such cases, it is sufficient in the above construction to store the previous value of this message field. We can sometimes decrease the message field by taking several steps of $`f_2`$ to simulate a single step of $`f_1`$. ∎

In case of one-dimensional systems, the “marching soldiers” scheme has the following strengthening of the original property saying that $`\zeta (x,t)`$ is independent of the order of updating. as in Example 1.1

###### Theorem 3.3.

Let $`𝔸_1=(,\mathrm{SS}_1,N,f_1)`$ be an arbitrary one-dimensional cellular automaton defined, as in Example 1.1, via a transition function $`g`$. Let the automaton $`𝔸_2=(,\mathrm{SS}_1\times R,N,f_2)`$ be defined as in the proof of Theorem 3.1. Let $`\eta `$ be an arbitrary asynchronous trajectory of $`𝔸_2`$. Let us define the functions $`\delta (x)`$, $`\overline{\eta }(x,u)`$ by $`\delta (0)=0`$, and

$`\delta (x+1)`$ $`=\delta (x)+\eta (x+1,0).\text{Age}\eta (x,0).\text{Age},`$
$`\overline{\tau }(x,t)`$ $`=\tau (x,t)+\delta (x),`$
$`\overline{\eta }(x,u)`$ $`=\zeta (x,u\delta (x)).\text{Cur}`$
for all $`u`$ of the form $`\overline{\tau }(x,t)`$. Also, let $`\overline{\eta }(x,\delta (x)1)=\eta (x,0).\text{Prev}`$. Then $`\tau (x,t)>0`$ implies with $`u=\overline{\tau }(x,t)1`$ that

$$\overline{\eta }(x,u+1)=g(\overline{\eta }(x1,u),\overline{\eta }(x,u),\overline{\eta }(x+1,u))$$

and all terms in this equation are defined.

The proof is straightforward verification. The theorem essentially says that from each asynchronous trajectory $`\eta `$ of $`𝔸_2`$, some synchronous trajectory $`\overline{\eta }`$ of $`𝔸_1`$ can be reconstructed as $`\overline{\eta }(x,u)=\zeta (x,u\delta (x)).\text{Cur}`$. The function $`\delta (x)`$ shows how much “ahead” or “behind” we are in simulating this trajectory when we start in $`\eta `$.

###### Remark 3.4.

This theorem fails in other neighborhood structures, namely in networks containing cycles: there, only certain initial configurations $`\eta (,0)`$ allow the construction of $`\delta (x)`$. In the ones that do not allow it, there is some inconsistency in the timing function $`\eta (x,0).\text{Age}`$ (a loop along which the sum of local increments of Age is not 0). In a connected network, this loop will imply that each cell can have only finitely many state changes, even in an infinite trajectory. ∎

## 4. Undecidability

###### Lemma 4.1.

Let us be given a one-dimensional commutative cellular automaton over the set of natural numbers, with “free boundary condition”, by a set of states $`\mathrm{SS}=\{0,\mathrm{},n1\}`$, transition functions $`g:\mathrm{SS}^3\mathrm{SS}`$ and $`g_0:\mathrm{SS}^2\mathrm{SS}`$ as in Example 1.1, with $`g(0,0,0)=0`$, $`g_0(1,s)=1`$ (for all $`s`$). The following problem is undecidable, as a function of $`n,g,g_0`$: Is there any synchronous trajectory of this cellular automaton, with $`\eta (x,0)=0`$ for all $`x`$ and $`\eta (0,t)=1`$ for some $`t>0`$?

###### Proof.

There is a standard construction to simulate Turing machines with such cellular automata, so the question reduces to the question whether an arbitrary Turing machine will halt when started on an empty tape. ∎

###### Lemma 4.2.

Let us be given a one-dimensional commutative cellular automaton over the set of natural numbers, with “free boundary condition”, by a set of states $`\mathrm{SS}=\{0,\mathrm{},n1\}`$, transition functions $`g:\mathrm{SS}^3\mathrm{SS}`$ and $`g_0:\mathrm{SS}^2\mathrm{SS}`$ as in Example 1.1.

The following problem is undecidable, as a function of $`n,g,g_0`$: Is there any trajectory of this cellular automaton, with $`\eta (0,0)=0`$ and $`\eta (0,t)=1`$ for some $`t>0`$?

Of course, once the automaton is commutative it does not matter whether the trajectory asked for is synchronous or asynchronous.

###### Proof.

From now on, without danger of confusion, let us write $`g(r,s)=g_0(r,s)`$ and forget about $`g_0`$. Let us be given a cellular automaton $`g`$ like in Lemma 4.1, with state set $`\mathrm{SS}=\{0,\mathrm{},n1\}`$. We construct a new cellular automaton over the set of states $`\mathrm{SS}^{}=\mathrm{SS}\{n\}`$, with the following transition function $`g^{}`$. Over states $`s<n`$, the functions $`g^{}`$ behave as $`g`$. Further, we have the following rules for $`g^{}`$ when at least one of the arguments is $`n`$.

$`(n,s)`$ $`g(0,0),`$
$`(s,n)`$ $`g(s,0)`$ $`\text{for }s<n,`$
$`(n,r,s)`$ $`n,`$
$`(r,n,s)`$ $`g(r,0,0)`$ $`\text{for }r<n,`$
$`(r,s,n)`$ $`g(r,s,0)`$ $`\text{for }r,s<n,`$
and $`(r,s,n)s`$, $`(r,s)r`$ in all remaining cases. By these rules, the symbol $`n`$ “sweeps” right and in its wake the rule $`g`$ will operate as if it had started from the a configuration of all 0’s. Thus, let $`\eta `$ be the synchronous trajectory of $`g`$ with $`\eta (x,0)=0`$ for all $`x`$. Then clearly if $`\eta ^{}`$ is any synchronous trajectory of $`g^{}`$ with $`\eta ^{}(0,0)=n`$ then for all $`t>0`$, for all $`xt`$ we have $`\eta ^{}(x,t)=\eta (x,t)`$.

Let us now apply the construction of the proof of Theorem 3.1 to $`g^{}`$ to obtain commutative rule $`g^{\prime \prime }`$ over the set of states $`\mathrm{SS}^{\prime \prime }=(\mathrm{SS}^{})^2\times \{0,1,2\}`$. We will prove that $`g^{\prime \prime }`$ has an asynchronous trajectory $`\eta ^{\prime \prime }`$ with $`\eta ^{\prime \prime }(0,0)=(n,0,0)`$ and $`\eta ^{\prime \prime }(0,u)=(1,1,0)`$ for some $`u`$, if and only if $`g`$ has a synchronous trajectory $`\eta `$ with $`\eta (0,x)=0`$ for all $`x`$ and $`\eta (0,u)=1`$ for some $`u`$. Since we know that the question whether this happens is undecidable from $`g`$, we will have proved that the question whether some cellular automaton has an asynchronous trajectory $`\eta `$ with $`\eta (0,0)=s_0`$ and $`\eta (0,u)=s_1`$ for some $`s_0s_1`$ is undecidable; this will complete the proof.

The “if” part: Suppose first that $`g`$ has a synchrounous trajectory $`\eta `$ with $`\eta (0,x)=0`$ for all $`x`$, and and $`\eta (0,u)=1`$ for some $`u`$. As mentioned above, then the synchronous trajectory $`\eta ^{}`$ of $`g^{}`$ has $`\eta ^{}(x,t)=\eta (x,t)`$ for all $`xt`$. Consider the synchronous trajectory $`\eta ^{\prime \prime }`$ of $`g^{\prime \prime }`$ started from $`\eta ^{\prime \prime }(x,0)=(n,0,0)`$ for all $`x`$. Then for all $`t>0`$ and all $`xt`$ we have

$$\eta ^{\prime \prime }(x,t)=(\eta ^{}(x,t),\eta ^{}(x,t1),tmod3)=(\eta (x,t),\eta (x,t1),tmod3).$$

Let $`v`$ be the first number $`>u+1`$ divisible by 3. We have

$$\eta ^{\prime \prime }(0,v)=(\eta (0,v),\eta (0,v1),0)=(1,1,0).$$

The “only if” part: Assume that $`\eta ^{\prime \prime }`$ is an asynchronous trajectory of $`g^{\prime \prime }`$ with $`\eta ^{\prime \prime }(0,0)=(n,0,0)`$ and $`\eta ^{\prime \prime }(0,w)=(1,1,0)`$ for some $`w`$. Then $`\tau ^{\prime \prime }(0,w)>0`$ and defining $`u=\overline{\tau }^{\prime \prime }(0,w)1`$, Theorem 3.3 implies

$$\overline{\eta }^{\prime \prime }(0,u+1)=g^{}(\overline{\eta }^{\prime \prime }(0,u),\overline{\eta }^{\prime \prime }(1,u)).$$

Applying theorem repeatedly, we obtain

$$\overline{\eta }^{\prime \prime }(0,v+1)=g^{}(\overline{\eta }^{\prime \prime }(x1,v),\overline{\eta }^{\prime \prime }(x,v),\overline{\eta }^{\prime \prime }(x+1,v))$$

or, if $`x=0`$, the same relation with the first argument of $`g^{}`$ omitted, for $`v=0,\mathrm{},u`$ and $`x\mathrm{min}\{v,(uv)\}`$. Now, if $`\eta ^{\prime \prime }(0,w)=(1,1,0)`$ then $`\overline{\eta }^{\prime \prime }(0,u+1)=1`$ while $`\overline{\eta }^{\prime \prime }(0,0)=n`$. We have just found that $`\overline{\eta }^{\prime \prime }(x,v)`$ develops according to $`g^{}`$ for $`v=0,\mathrm{},u`$ and $`x\mathrm{min}\{v,(uv)\}`$. As discussed above, therefore $`\overline{\eta }^{\prime \prime }(0,u+1)=1`$ if and only if $`g`$ computes 1 at $`(0,u+1)`$ from an all-0 initial configuration. ∎

###### Proof of Theorem 1.3.

Let the local state space be the set of integers $`\mathrm{SS}=\{0,\mathrm{},n+2\}`$. Let $`g:\mathrm{SS}_0^3\mathrm{SS}_0`$ and $`g_0:\mathrm{SS}_0^2\mathrm{SS}_0`$ be the rules for a commutative cellular automaton transition rule with state set $`\mathrm{SS}_0=\{0,\mathrm{},n1\}`$. We define the transition function $`f`$. We will write $`f(x,y,z)=y^{}`$ as $`(x,y,z)y^{}`$. We require

(4.1) $`(s,n,0)`$ $`n+1,`$
(4.2) $`(s,n,1)`$ $`n+2,`$
(4.3) $`(r,s,t)`$ $`g_0(s,t),`$ $`\text{for all }rn,r,s<n,r1,`$
(4.4) $`(r,s,t)`$ $`g(r,s,t)`$ $`\text{for all }r,s,t<n,`$
(4.5) $`(r,s,t)`$ $`g(r,s,0)`$ $`\text{for all }r,s<n,tn,`$
and $`(r,s,t)s`$ in all remaining cases. Let us show that $`f`$ has invariant histories if and only if $`g`$ has no asynchronous trajectory $`\eta _0`$ over $`=_+`$ with $`\eta _0(0,0)=0`$ and $`\eta _0(0,t)=1`$ for some $`t`$. Assume first that $`g`$ has such a trajectory. Let us define the initial configuration $`\xi `$ of $`f`$ as $`\xi (x)=n`$ if $`x=1`$ and 0 otherwise. We may apply rule (4.1) first to get $`\eta (1,1)=n+1`$. Or, we may apply rules (4.3),(4.4),(4.5) first to cells $`x>0`$ on the right repeatedly. Sooner or later we have $`\eta (0,t)=1`$, which allows $`\eta (1,t+1)=n+2`$ by rule (4.2) in the next step. Thus, depending on the order of rule application, we obtained in cell $`1`$ the sequence $`n,n+1`$ or $`n,n+2`$.

Suppose now that $`g`$ has no such trajectory and let $`\xi `$ be an arbitrary configuration of $`f`$. Each occurrence of a state $`n`$ remains such an occurrence. On segments between them, the commutative rule $`g`$ works. The only other transitions possible are $`(r,n,0)n+1`$ and $`(r,n,1)n+2`$. Assume $`\eta (x,0)=n`$ and consider the sequence of different values in $`\eta (x+1,t)`$. Let us show that 0 and 1 cannot both occur in this sequence and hence only one of the transitions is possible. Indeed, if 0 occurs before 1 then our assumption about $`g`$ excludes the occurrence of 1 in the sequence any later. If 1 occurs in the sequence before 0 then our rules (in particular, (4.3)) do not allow any change of the state of $`x+1`$ after that. ∎

## Acknowledgments

I thank Robert Solovay for pointing out several errors in the first version, Wayne Snyder for calling my attention to  and the anonymous referee for his careful reading and many corrections.


format: text

text:
# Casimir and pseudo-Casimir interactions in confined polyelectrolytes

## 1 Introduction

Charged polymers have received much attention over the years for their indisputable importance in industrial applications as well as for their role in fundamental colloid science (for comprehensive introductions see ). Apart from the study of polymers in the bulk, their interactions with surfaces, as well as polymer-mediated interactions between surfaces have been studied on different levels . Though polyelectrolyte interaction with charged surfaces is quite similar to the general interaction of polymers with neutral surfaces , the long range nature of the Coulomb potential introduces additional features and difficulties that make these problems even more difficult to handle.

Most of the work on inhomogeneous polyelectrolytes was done on the mean-field level and the effect of thermal fluctuations has not been considered at all. The mean-field picture of the behavior of a polyelectrolyte chain confined between two oppositely charged surfaces that emerged from previous studies  (see Fig. 1) lets one think that at small intersurface separations the polyelectrolyte mediated interations are quite similar to those in the case of neutral polymers and decay algebraically with separation . At a critical value of the intersurface spacing there is a transition of the confined polyelectrolyte chain from a monomodal to a bimodal monomer density distribution. This means that at sufficiently large intersurface separations the polymer chain is on the average partly electrostatically adsorbed to both surfaces and the remaining part of the chain bridges the region between the surfaces. The consequence of these bridging interactions is an exponentially decaying attraction between the surfaces that persists for all separations greater than the critical value . In this respect the charged chain between oppositely charged surfaces behaves on the mean-field level a lot like a neutral chain with a finite adsorption energy to the two bounding surfaces.

Recent theoretical studies  however point to conclusion that thermal fluctuation effects in Coulombic systems, ignored on the mean-field level, can hardly ever be considered as negligible. In the limit of small electrostatic coupling it has been known for a while that thermal fluctuations in an inhomogeneous Coulomb fluid lead to zero order van der Waals (or equivalently Casimir) forces . More recent studies of thermal effects in soft matter in general  are very much consistent with these findings. It has been established that thermally driven fluctuations of order parameters of confined systems also lead to long range pseudo-Casimir forces , quite similar in nature to ordinary Casimir forces except that instead of being due to electromagnetic field fluctuations they are a consequence of a general order parameter fluctuations (thus pseudo-Casimir interactions). In view of this we suspect that similar phenomena should also show up in the study of fluctuation effects of confined inhomogeneous polyelectrolytes.

In order to prove this conjecture we write the partition function for a polyelectrolyte chain confined between two oppositely charged planar surfaces in such a way that we were able to expand it to quadratic order in the relevant order parameters, in this case the local electric field and the replicated polymer density field (see below), around the mean-field profile and evaluate on this level the effect of thermal fluctuations. We note that thermal fluctuations decouple into two separate contributions: the Casimir contribution due to electric field fluctuations and the pseudo-Casimir contribution due to the replicated polymer density field fluctuations. We show that the effect of thermal fluctuations depends on only one coupling parameter $`\alpha _{MF}=\frac{4\pi }{3}\frac{\sigma \mathrm{}}{\tau }`$, where $`\sigma `$ is the surface charge density, $`\tau `$ is the line charge density of the polymer chain and $`\mathrm{}`$ is the Kuhn’s segment length. Depending on the value of this parameter the behavior of the system can be successfully described on the mean-field level ($`\alpha _{MF}1`$) or has to be ammended by taking into account the fluctuation terms ($`\alpha _{MF}1`$). The most important consequence of the fluctuation effects seems to be the long range pseudo-Casimir force, with an anomalous algebraic dependence on the intersurface separations with an exponent $`1.7`$, due to the confined fluctuations of the polymer density field that decays slower than the usual Casimir (or equivalently van der Waals) force.

The plan of the paper is as follows. In order to investigate the fluctuation effects we first write the free energy of a confined charged polymer chain between two surfaces in the form of an $`𝒪(n)`$ scalar field theory. From here the free energy is obtained via the standard $`n0`$ replica trick . Since the mean-field is given by the saddle-point, we then expand this field theory to the second order around the saddle-point in both the electrostatic field and the replicated polymer density field. Calculation of the quadratic corrections to the mean-field free energy that can be performed exactly gives the Casimir and pseudo-Casimir interactions mediated by the confined electrostatic and confined replicated polymer density fields. We investigate the dependence of the regularized fluctuation interaction free energy on the intersurface separation in the limit of small and large spacings compared to the critical separation.

## 2 Mean-field theory

As a point of departure we shall take a model confined polyelectrolyte system, composed of two surfaces with a specified surface charge density ($`\sigma `$) and a confined, oppositely charged single polyelectrolyte chain of length $`N\mathrm{}`$, described with an Edwards hamiltonian

$$[𝐫(s)]=\frac{3}{2\mathrm{}}_0^N\mathrm{}\dot{𝐫}^2(s)𝑑s+\frac{1}{2}_0^N\mathrm{}_0^N\mathrm{}V(𝐫(s),𝐫(s^{}))𝑑s𝑑s^{}.$$
(1)
We assume that the intersegment interactions are given by the unscreened Coulombic form $`V(𝐫,𝐫^{})=\tau ^2/4\pi ϵϵ_0|𝐫𝐫^{}|`$, where $`\tau `$ is the charge per unit length. This is a bit artificial since we completely ignore the presence of other mobile ionic species such as counterions and salt ions, but represents a clean, tractable model system where fluctuation effects can be studied in detail.

Performing now the Hubbard-Stratonovich transformation for the pair interaction potential, going to the grand canonical ensemble with a fixed chemical potential $`\mu `$ for the monomers and adequatly removing the closed loop polymer configurations from the partition function via the $`n0`$ limit of an $`𝒪(n)`$ scalar field theory with the field $`𝝍=(\psi _1,\psi _2\mathrm{}\psi _n)`$, we remain with the de Gennes - des Cloizeaux representation for the free energy of a charged polyelectrolyte chain

$$(\mu ,\beta )=\underset{n0}{lim}\frac{1}{n}\mathrm{log}𝒵(\mu ,\beta )=lim_{n0}\frac{1}{n}\mathrm{log}𝒟\varphi (𝐫)𝒟𝝍(𝐫)\mathrm{exp}\left(\beta d^3𝐫[\varphi (𝐫),\psi _i(𝐫)]\right)$$
(2)
where

$$\beta [\varphi (𝐫),\psi _i(𝐫)]=\frac{1}{2}\frac{\mathrm{}^2}{6}_i\left(\mathbf{}\psi _i\right)^2+\frac{1}{2}\mu _i\psi _i^2+\frac{1}{2}\beta ϵϵ_0\left(\mathbf{}\varphi \right)^2+i\beta \tau _i\varphi \psi _i^2.$$
(3)
In writing Eq. 2 we omitted the determinant $`\mathrm{det}\left(\beta ϵϵ_0^2\right)`$ from the denominator. As is usual on the one-loop expansion level, this term is exactly cancelled by the zero order van der Waals term that one should add to the final total interaction free energy. This point is thoroughly explained in . See also below.

The path from here will be to obtain the mean-field solution of the model defined by Eq. 2 via the saddle-point of the $`𝒪(n)`$ scalar field theory along the lines of  and then to evaluate the contribution of quadratic fluctuations around the saddle-point to the thermodynamic properties of the system.

Evaluation of the saddle-point from the field hamiltonian Eq. 3 is trivial. One notices immediatley that all the $`\psi _i`$ satisfy the same saddle-point equation and thus $`\psi _i\psi `$. With this one can write the dependence of the hamiltonian Eq. 3 on $`n`$ explicitly  and the limit $`n0`$ can be obtained straightforwardly. One remains with

$$\beta _0=\underset{n0}{lim}\frac{1}{n}\mathrm{log}𝒵_0=\frac{1}{2}d^3𝐫\left(\frac{\mathrm{}^2}{6}\left(\mathbf{}\psi \right)^2+\mu \psi ^2+\beta ϵϵ_0\left(\mathbf{}\varphi \right)^2+i\beta \tau \varphi \psi ^2\right),$$
(4)
where we introduced the mean polymer density field $`\psi `$ so that the local monomer density is $`\rho =\psi ^2`$ and the local electrostatic mean potential $`\phi =i\varphi `$. In the planparallel geometry considered here all the fields depend only on the transverse coordinate $`z`$, and the inhomogeneities are confined to this direction only $`\psi =\psi (z),\phi =\phi (z)`$. The saddle-point equations of the hamiltonian Eq. 3 can thus be written in the form

$`{\displaystyle \frac{\mathrm{}^2}{6}}{\displaystyle \frac{d^2\psi (z)}{dz^2}}=\mu \psi (z)+2\beta \tau \phi (z)\psi (z),ϵϵ_0{\displaystyle \frac{d^2\phi (z)}{dz^2}}=\tau \psi ^2(z).`$ (5)
These mean-field equations are exactly equivalent to those derived previously  if one takes into account that the polyelectrolyte chain is the only mobile charge in the system (no counterions and no salt). The first equation of Eq. 5 is the equation for the density field $`\psi (z)`$ of the polymer in an inhomogeneous external field $`\phi (z)`$ while the second one is a Poisson - Boltzmann equation for the mean electrostatic potential $`\phi (z)`$ of a charged polymer chain with charge density $`\rho (z)=\tau \psi ^2(z)`$.

The second equation of Eq. 5 can be solved explicitly and its solution can be manipulated  to yield a limiting form valid close to the boundaries at $`z\pm a`$ where the potential is largest

$$ϵϵ_0\phi (z)=\tau _a^{+a}|zz^{}|\psi ^2(z^{})𝑑z^{}\tau |z|+𝒪(z^2).$$
(6)
This expression can be obtained rather straightforwardly by Taylor expanding the solution around $`z=\pm a`$. In the opposite limit, $`z0`$, we can derive

$$ϵϵ_0\phi (z)\tau _a^{+a}|u|\psi ^2(u)𝑑u\tau \psi ^2(0)z^2+𝒪(z^4).$$

Since the mean-field is largest close to the boundaries we will use the approximate expression valid strictly only close to $`z=\pm a`$ in the whole interval. This approximation works extremely well everywhere  except close to the origin where $`z`$ is small anyway. The reason for introducing this additional but inessential approximation is that it also helps in transforming the first of Eqs. 5 from a non-linear Landau-Ginzburg type equation into a linear equation, which can be solved analytically and explicitly in terms of Airy’s functions .

The electroneutrality demands that at the bounding surfaces ($`z=\pm a`$) we have $`ϵϵ_0\frac{\phi }{𝐧}=\tau \frac{N}{S}\sigma `$, where $`𝐧`$ is the boundary surface normal. Because of the impenetrability of the boundaries to the polymer chain we should also have $`\psi (z=\pm a)=0`$. The unnormalized solution of the first mean-field equation Eq. 5 can now be obtained explicitly in the form

$$\psi (x)Ai(y_0x)Bi^{}(y_0)Bi(y_0x)Ai^{}(y_0),$$
(7)
where we introduced the dimensionless variables $`x=\lambda _B^{1/3}z`$, $`x_0=\lambda _B^{1/3}a`$, $`y_0=6\mu \lambda _B^{2/3}/\mathrm{}^2`$, $`\lambda _B=\frac{12\beta \tau ^2}{ϵϵ_0\mathrm{}^2}`$ and $`Ai(x)`$ and $`Bi(x)`$ are the standard Airy functions. The dependence $`y_0=y_0(x_0)`$ is obtained from the vanishing density field boundary conditions at the two bounding surfaces. We derive  that $`y_0(x_0)`$ grows linearly with $`x_0`$ for large $`x_0`$, and therefore $`u_0=y_0x_0`$ asymptotically approaches a constant $`u_0^{\mathrm{}}2.34`$, while for small $`x_0`$, $`y_0`$ is negative and behaves approximately like $`y_0(\pi /2x_0)^2`$.

The main consequence of the mean-field equations is that for very small intersurface separations the chain is desorbed with a monomodal density distribution . In this regime of intersurface separations the interactions between the surfaces are repulsive, decaying albebraically with plate separation. At a critical value of the intersurface separation $`x_0=1.986`$ the electrostatic attraction causes the chain to adsorb to both surfaces creating a polymer bridge leading in its turn to bridging attraction between the surfaces. This attraction, though exponentially small, persists for all separations between the surfaces larger than the critical. The intersurface pressure is thus repulsive for $`x_0<1.986`$ and attractive for $`x_0>1.986`$. For details see .

## 3 Fluctuations

We now turn to the contribution of fluctuations around the saddle point to the free energy. First we expand the hamiltonian Eq. 3 to the second order in deviations from the mean-fields Eq. 5. Our treatment of the fluctuation effect will thus be based on a quadratic (one-loop) expansion. The second functional derivatives of the field hamiltonian can be assembled into a Hessian of the fields of Eq. 3

$$=\left[\begin{array}{cc}\frac{\mathrm{}^2}{6}^2+\mu +2i\beta \tau \varphi & 2i\beta \tau \psi _n\\ 2i\beta \tau \psi _n& \beta ϵϵ_0^2\end{array}\right],$$
(8)
where $`\phi `$ and $`\psi _n`$ are the solutions of the mean-field equations Eq. 5. While integrating out the mean replicated polymer density field and the mean electrostatic field fluctuations we notice that the functional integral over $`\psi _i`$ fluctuations can be evaluated directly, giving an explicit dependence of the fluctuational partition function on $`n`$. This is similar to the general case treated by Emery . The limit $`n0`$ can now be evaluated explicitly leading to the following contribution of the fluctuations around the mean-fields in the quadratic order to the free energy

$`\beta _2`$ $`=`$ $`\underset{n0}{lim}\frac{1}{n}\mathrm{log}𝒵_2=\frac{1}{2}\mathrm{lndet}\left(^2\right)+\frac{1}{2}\mathrm{lndet}\left(^2+\frac{6}{\mathrm{}^2}\mu +12\frac{\beta \tau }{\mathrm{}^2}\phi (𝐫)\right)+`$ (9)
$`+`$ $`12{\displaystyle \frac{(\beta \tau ^2)}{ϵϵ_0\mathrm{}^2}}{\displaystyle d^3𝐫d^3𝐫^{}𝒢(𝐫,𝐫^{})𝒢_0(𝐫,𝐫^{})\psi (𝐫)\psi (𝐫^{})},`$
where $`\phi (𝐫)`$ and $`\psi (𝐫)`$ are the solutions of the mean-field equations Eq. 5. The two Green’s functions $`𝒢(𝐫,𝐫^{})`$ and $`𝒢_0(𝐫,𝐫^{})`$ are the functional inverses defined as

$`\left(^2+{\displaystyle \frac{6}{\mathrm{}^2}}\mu +12{\displaystyle \frac{\beta \tau }{\mathrm{}^2}}\phi (𝐫)\right)𝒢(𝐫,𝐫^{})`$ $`=`$ $`\delta ^3(𝐫𝐫^{})`$
$`\left(^2\right)𝒢_0(𝐫,𝐫^{})`$ $`=`$ $`\delta ^3(𝐫𝐫^{}).`$ (10)
The first term in Eq. 9 comes from the confined electric field fluctuations and is in fact the standard zero-order Casimir interaction energy  (see below). The second and the third terms are due to the confined polymer density field fluctuations and thus correspond to the pseudo-Casimir interactions. The free energy Eq. 9 is obtained already after the cancellation referred to in the discussion following Eq. 3 and thus represents the total interaction free energy including the zero order van der Waals contribution. The latter is evaluated for a slab of thickness $`2a`$ with boundaries impenetrable to electrostatic fields.

Let us first deal with the two functional determinants in Eq. 9. Since the mean-field solution depends only on the transverse coordinate $`z`$ we need to investigate the functional determinant of an operator of the type $`\frac{d^2}{dz^2}+V(z)`$ that can be evaluated in two different ways. One can first of all find the eigenvalues $`\lambda _n`$ of this operator with the boundary condition $`f_n(z=\pm a)=0`$. The functional determinant comes out as the product of these eigenvalues $`\mathrm{\Pi }_n\lambda _n`$. The other approach would take into account the van Vleck identity  involving $`f_0(z)`$, the eigenfunction with zero eigenvalue and with boundary conditions $`f_0(z=a)=0`$ and $`f_0^{}(z=a)=1`$. The logarithm of the functional determinant can thus be written equivalently (up to an irrelevant - in this context - additive constant) in two different ways

$$\mathrm{ln}\mathrm{det}\left(\frac{d^2}{dz^2}+V(z)\right)=\underset{n}{}\mathrm{ln}\lambda _n=\mathrm{ln}\left[f_0(z=a)\right].$$
(11)
In what follows we will use both representations together with the Fourier decomposition in the $`(x,y)`$ plane and will find both equally inadequate when approaching the critical intersurface separation $`x_0=1.986`$, either from above or below, where the mean polymer density field goes through a monomodal to bimodal transition.

Of the first two terms in Eq. 9 the first one is trivial to evaluate either via the van Vleck relation or the eigenvalue product representation. First of all we introduce the Fourier decomposition in the $`(x,y)`$ plane with a wave vector $`𝐐`$, since we only need to take into account the boundaries in the $`z`$ direction. This leads to the standard zero-order (thermal) Casimir interaction in the form

$$\frac{1}{2}k_BT_𝐐\mathrm{ln}\mathrm{det}\left(\frac{d^2}{dz^2}+Q^2\right)=\frac{Sk_BT}{4\pi }_0^{\mathrm{}}Q𝑑Q\mathrm{ln}\left(\frac{\mathrm{cosh}Qa\mathrm{sinh}Qa}{Q}\right),$$
(12)
where $`S`$ is the total area of the interacting boundary surfaces.

The second term is trickier. First of all we again introduce the Fourier decomposition in the $`(x,y)`$ plane with a wave vector $`𝐐`$ and invoke the approximate solution Eq. 6. Using again the van Vleck identity we derive

$$\frac{1}{2}k_BT_𝐐\mathrm{ln}\mathrm{det}\left(\frac{d^2}{dz^2}+Q^2+\frac{6}{\mathrm{}^2}\mu +12\frac{\beta \tau }{\mathrm{}^2}\phi (z)\right)=\frac{Sk_BT}{4\pi }_0^{\mathrm{}}Q𝑑Q\mathrm{ln}\left(C(Q,a)S(Q,a)\right),$$
(13)
where we decomposed the determinant of the operator into two separate contributions stemming from symmetric and antisymmetric modes

$`C(Q,z)`$ $`=`$ $`Bi^{}(y(Q))Ai(y(Q)\lambda _B^{1/3}z)Ai^{}(y(Q))Bi(y(Q)\lambda _B^{1/3}z)`$
$`S(Q,z)`$ $`=`$ $`Bi(y(Q))Ai(y(Q)\lambda _B^{1/3}z)Ai(y(Q))Bi(y(Q)\lambda _B^{1/3}z),`$ (14)
where $`y(Q)=\lambda _B^{2/3}\left(Q^2+\frac{6\mu }{\mathrm{}^2}\right)`$ while $`Ai(x)`$ and $`Bi(x)`$ are the standard Airy functions . The Airy function solutions are of course a consequence of the approximation Eq. 6 for the mean electrostatic potential. The free energy Eq. 13 still needs to be properly regularized by subtracting the infinite bulk and surface terms, see below.

In what follows we introduced also these dimensionless variables: $`q^2=\lambda _B^{2/3}Q^2`$ and $`u(q,x,x_0)=y(q,x_0)x=q^2+y_0(x_0)x=u_0(x_0)+q^2+x_0x`$. Here $`y_0(x_0)`$ is the largest zero of the symmetric function $`C(y_0(x_0),x_0)`$ which corresponds exactly to the zero of the mean-field solution Eq. 7 at the two bounding surfaces .

The third term in Eq. 9 appears to be the most difficult one to evaluate but fortunately it is possible to prove that within the approximation Eq. 6 its contribution to the fluctuation free energy $`_2`$ is in fact negligible compared to the other two, i.e. Eqs. 13 and 14. Taking into account the approximate form for the mean electrostatic potential Eq. 6 the last term in Eq. 9 can be evaluated explicitly leading to

$$12\frac{(\beta \tau ^2)}{ϵϵ_0\mathrm{}^2}d^3𝐫d^3𝐫^{}𝒢(𝐫,𝐫^{})𝒢_0(𝐫,𝐫^{})\psi (𝐫)\psi (𝐫^{})=\frac{\lambda _B}{(ϵϵ_0)}_a^{+a}𝑑z|z|𝒢(z,z^{},Q=0),$$
(15)
where we took into account that the mean-field solution $`\psi (𝐫)`$ depends only on the transverse coordinate $`z`$ and that again we can introduce the Fourier decomposition in the $`(x,y)`$ plane with a wave vector $`𝐐`$. The magnitude of the coupling contribution Eq. 15 turns out to be much smaller than the other two, Eqs. 12, 13, and can thus be safely ignored in numerical computations.

The fluctuation free energy, $`_2`$ which is thus composed only of the contributions from the fluctuational determinants of the electrostatic field and the replicated polymer density field, still contains the divergent bulk and surface parts together with the interaction terms . $`_2`$ is thus formally infinite. This infinity can be regularized by subtracting the infinite bulk and surface contributions  and thus in forming the regularized fluctuation free energy

$$𝒲_2(a)=_2(a)_2^{\mathrm{}}(a).$$

Here $`_2^{\mathrm{}}(a)`$ stands for the form of the fluctuation free energy for large values of the argument but evaluated at a finite $`a`$ .

Using now the standard form for the mean-field pressure  and adding to it the fluctuation contribution, one can write for the total pressure between the boundary surfaces

$`p(a)`$ $`=`$ $`p_{MF}+p_{FL}[\delta \varphi ]+p_{FL}[\delta \psi ]=`$ (16)
$`=`$ $`\frac{1}{2}k_BT\mu \psi ^2(z=0)\frac{k_BT}{8\pi }\frac{\zeta (3)}{(2a)^3}+p_{FL}(a),`$
where the first term on the r.h.s. derives from the mean-field solution, the second one derives from the Casimir term (stemming from the electrostatic field fluctuations, first term in Eq. 9 and the third one derives from the pseudo-Casimir terms (stemming from the replicated density field fluctuations, second term in Eq. 9 with

$$p_{FL}(a)=\frac{1}{S}\left(\frac{𝒲_2(a)}{(2a)}\right).$$
(17)
The next section deals with the evaluation and regularization of the fluctuation terms in Eq. 16.

## 4 Results - large intersurface separation

We will first have to regularize the divergent expression Eq. 13 by subtracting the divergent terms explicitly in the limit of large $`x_0`$ (to be quantified below). When $`x_0`$ is large enough the argument $`y=u+x_0`$ of the Airy functions in Eq. 14 becomes large too and we can use the asymptotic expansion formulas  for the Airy functions to explicitly find the divergent terms in $`_2^{\mathrm{}}(a)`$. The regularized free energy is then

$`𝒲_2(x_0)={\displaystyle \frac{Sk_BT}{4\pi }}{\displaystyle _{u_0}^{\mathrm{}}}𝑑u`$
$`\mathrm{ln}\left({\displaystyle \frac{\left(Bi^{}(u+x_0)Ai(u)Ai^{}(u+x_0)Bi(u)\right)\left(Bi(u+x_0)Ai(u)Ai(u+x_0)Bi(u)\right)}{\left(Bi_L^{}(u+x_0)Ai(u)Ai_L^{}(u+x_0)Bi(u)\right)\left(Bi_L(u+x_0)Ai(u)Ai_L(u+x_0)Bi(u)\right)}}\right),`$
(18)
where $`Ai_L`$, $`Bi_L`$, $`Ai_L^{}`$ and $`Bi_L^{}`$ are the lowest order asymptotic expressions for the Airy functions . The next order in the asymptotic expansion of the Airy function is $`Ai_{L2}(t)0.07(\frac{2}{3}t^{3/2})^1Ai_L(t)`$ and in order for it to be negligible, $`t`$ has to be large enough, typically larger then $`1.5`$. This leads to the condition $`u_0(x_0)+x_0>1.5`$, which yields roughly that $`x_0`$ has to be larger than $`3`$ if our calculation is to be accurate to within a few percents. For smaller values of the intersurface separation the form of the regularized interaction free energy Eq. 18 does not work.

Fig. 2 shows the dependence of the dimensionless fluctuational free energy $`𝒲_2`$ on the dimensionless separation $`x_0`$. As already stated the regularization procedure Eq. 18 can only be used for sufficiently large intersurface spacing. The log-log plot (see the inset of Fig. 2) reveals that asymptotically the fluctuational contribution with an algebraic decay prevails. From the inset of Fig. 2 we read that the free energy decays asymptotically as $`x_0^{1.7}`$. The algebraic decay of the free energy and concommitantly of the pressure, indicates that the fluctuation induced force has a much longer range than the mean field force which decays approximately exponentially . Also it has a longer range than the zero-order Casimir (or equivalently the zero order van der Waals) interaction free energy that decays as $`x_0^2`$.

There is also no obvious connection between this result and the behavior of the fluctuation part of the interaction free energy of a system composed of unconnected counterions treated on the same level of approximmation . In that case the fluctuation free energy decays as $`x_0^2\mathrm{ln}x_0`$, which is much faster than $`x_0^{1.7}`$.

## 5 Results - small intersurface separation

For small plate separations the above procedure cannot work because in that case we can not write down explicitly the asymptotic form of the Airy functions that would be valid in the whole range of integration in Eq. 13. We found out that the best way to approach the evaluation of the fluctuation determinant in this limit is via a perturbation analysis using the eigenvalue representation of the functional determinant Eq. 11.

The perturbative treatment of the evaluation of Eq. 13 can be introduced as follows. The differential equation we have to solve is

$$\left(\frac{d^2}{dx^2}+q^2+y_0(x_0)+|xx_0|\right)\psi =0.$$
(19)
The last term can be treated as a linear perturbation potential $`V(x)=|xx_0|`$. What we will do now is to evaluate the determinant of the operator in Eq. 19 in the form of a perturbation expansion in $`V(x)`$ up to and including the first order. The order zero of the perturbation expansion of the determinant of the operator Eq. 19 gives

$$\mathrm{det}\left(\frac{d^2}{dx^2}+q^2+y_0\right)=\psi _0^0(x_0)=\frac{\mathrm{sinh}2\sqrt{q^2+y_0}x_0}{2\sqrt{q^2+y_0}}.$$
(20)
The final integration over the wavevector $`𝐐`$ has to be divided into two parts because $`y_0`$ is negative, leading to

$$𝒲_2^0=\frac{Sk_BT\lambda _B^{2/3}}{16\pi x_0^2}\left[_0^\pi w𝑑w\mathrm{ln}\frac{\mathrm{sin}w}{w}\frac{\zeta (3)}{4}\right].$$
(21)
To obtain the first order term in the perturbation expansion we shall use the product of the eigenvalues representation of the determinant of the operator which leads to the following first order form of the free energy

$$_2^1=\frac{k_BT\lambda _B^{2/3}S}{4\pi }_0^{\mathrm{}}q𝑑q\underset{N}{}\mathrm{ln}\lambda _N^1$$
(22)
where $`\lambda _N^1`$ represents the first order correction to the eigenvalue of the differential equation Eq. 19. This correction can be evaluated by standard formulas of the perturbation theory  leading to the following result

$$\lambda _N^1=\lambda _N^0\frac{<N|V(x)|N>}{<N|N>}=\lambda _N^0\frac{_0^{2x_0}|xx_0|\mathrm{sin}^2\left(\frac{N\pi }{2x_0}x\right)𝑑x}{_0^{2x_0}\mathrm{sin}^2\left(\frac{N\pi }{2x_0}x\right)𝑑x}=\lambda _N^0+\frac{x_0}{2},$$
(23)
so the sum of the first order corrections to the log of the determinant of the operator Eq. 22 is

$$\underset{N}{}\mathrm{ln}\lambda _N^1\underset{N}{}\mathrm{ln}\lambda _N^0+\frac{x_0^3}{\pi ^2}\underset{N}{}\frac{1}{N^21+(\frac{2x_0}{\pi })^2(q^2+\frac{x_0}{2})}.$$
(24)
Since the first term in this expression has already been calculated, we now focus on the second one. This term is obviously divergent. If we calculated the bulk contribution, there would be an integral over $`N`$ instead of the discrete sum, which is a consequence of the boundary conditions. The regularized free energy is thus  obtained as the difference between the sum and the integral of the same expression. After some algebra the result of this manipulation is

$`𝒲_2^1={\displaystyle \frac{Sk_BT\lambda _B^{2/3}x_0}{16\pi }}[`$ $`\frac{1}{2}`$ $`{\displaystyle _0^{\sqrt{1\frac{2x_0^3}{\pi ^2}}}}𝑑w\left(\pi \left(\mathrm{coth}\pi w\mathrm{cot}\pi w1\right)\mathrm{ln}{\displaystyle \frac{1+w}{1w}}+2\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{t}\mathrm{a}\mathrm{n}{\displaystyle \frac{1}{w}}\right)+`$
$`+`$ $`\frac{1}{2}`$ $`{\displaystyle _{\sqrt{1\frac{2x_0^3}{\pi ^2}}}^{\mathrm{}}}dw(\pi (\mathrm{coth}\pi w1)+2(\mathrm{arctan}{\displaystyle \frac{1}{w}}{\displaystyle \frac{1}{w}}))].`$ (25)
Even without really calculating the second order we can nevertheless say something about its magnitude. The maximal difference between the first and the second order perturbation terms can be obtained in standard notation as

$$\delta \lambda _N^1\frac{1}{\mathrm{\Delta }\lambda }\left[<N|V(x)^2|N>\left(<N|V(x)|N>\right)^2\right],$$
(26)
where $`\mathrm{\Delta }\lambda `$ is the minimal difference between $`\lambda _N`$ and any other $`\lambda _M`$. A short calculation then yields $`\frac{\delta \lambda ^1}{\lambda ^1}\frac{x_0^3}{14.8}`$. As expected, the reliability of the first order result drops as we move towards larger $`x_0`$ i.e. towards the mean field phase transition from bellow.

Fig. 3 shows the first and the second perturbation order of the regularized fluctuation free energy at small intersurface separations. The pressure obviously changes sign for sufficiently large separation ($`x_0`$ larger than $`1.68`$ on the Fig. 3). For small enough separations the free energy scales approximately as $`x_0^2`$ which is the same scaling form as for the mean-field part. It seems to us that the change in sign of the pressure close to the transition point $`x_0=1.986`$ is not due entirely to artifacts of the first order perturbation theory for sufficiently large intersurface separations, but could concievably survive all higher order perturbation terms. Unfortunately all approximation schemes that we could think of broke down close to the transition point and we thus have to leave the important question about the behavior of the fluctuation component of the pressure close to the transition unanswered.

## 6 Discussion

We have described all the contributions to the fluctuational free energy of the system. To evaluate the fluctuational part of the pressure between the two surfaces the derivatives of the regularized free energy with respect to $`x_0`$ have been performed numerically. The total prerssure can thus be obtained by simply taking the sum of mean-field and fluctuational contributions Eq. 16. Fig. 4 shows the dimensionless total pressure as a function of the dimensionless separation $`x_0`$ which can be written, Eq. 16, as

$`p(x_0)`$ $`=`$ $`p_{MF}+p_{FL}[\delta \varphi ]+p_{FL}[\delta \psi ]=`$ (27)
$`=`$ $`{\displaystyle \frac{k_BT}{2}}{\displaystyle \frac{\lambda _B\mathrm{}^2\sigma }{8\tau }}y_0(x_0)\stackrel{~}{\psi }^2(0){\displaystyle \frac{k_BT\zeta (3)\lambda _B}{64\pi x_0^3}}\lambda _B^{1/3}{\displaystyle \frac{d}{d(2x_0)}}\left({\displaystyle \frac{𝒲_2(x_0)}{S}}\right)=`$
$`=`$ $`p^{}\left[y_0(x_0)\stackrel{~}{\psi }^2(0){\displaystyle \frac{\zeta (3)}{4\alpha _{MF}x_0^3}}+{\displaystyle \frac{f(x_0)}{\alpha _{MF}}}\right]`$
where $`f(x_0)`$ is a function of $`x_0`$ only , $`\stackrel{~}{\psi }(0)`$ is the normalized dimensionless $`\psi `$, $`\alpha _{MF}=\frac{4\pi }{3}\frac{\sigma \mathrm{}}{\tau }`$ and $`p^{}=\frac{k_BT\lambda _B\alpha _{MF}}{16\pi }=\frac{\sigma ^2}{2ϵϵ_0}`$. The parameter $`\alpha _{MF}`$ depends on the charge density on the plates $`\sigma `$, the linear charge density on the polyelectrolyte chain $`\tau `$ and on the segment length of the polyelectrolyte $`\mathrm{}`$. If we put in experimentally reasonable values of $`\sigma e_0/1,35nm^2`$ then $`\alpha _{MF}`$ is about 100 for DNA (a relatively stiff polyelectrolyte) and around 4 for hyaluronic acid (a relatively flexible polyelectrolyte). Obviously the larger $`\alpha _{MF}`$ the less important will the fluctuation effects be and the more accurate will be the mean-field result.

Fig. 4 shows the complete pressure Eq. 27 as a function of the dimensionless separation for different values of the coupling parameter $`\alpha _{MF}`$. As already stated we were unable to analyse the behavior of the interaction pressure close to the mean-field transition point. Our failure could indicate that close to the transition point our harmonic approach to the fluctuations breaks down. This would not be totally unexpected since we are dealing with a continuous transition in the distribution of the polyelectrolyte density that apparently has all the attributes of a second order phase transition, with large fluctuations close to the transition point. We will however not pursue this line of thought here.

It is clear from Fig. 4 that the effect of the fluctuations on the intersurface force can be profound. First of all the fluctuations modify the force at small separations making it less repulsive if compared to the mean-field case. Since both the mean-field as well as the fluctuation contribution to the free energy scale as $`x_0^2`$ (but with a different sign) it appears that the effect of the fluctuations is simply to renormalize the magnitude of the pressure but not its scaling form. If one would thus try to fit the complete $`p(a)`$ with $`p_{MF}(a)`$ one would have to introduce a smaller effective charge on the surfaces for the fit to make sense.

At large separations the behavior of $`p(a)`$ is completely different from the behavior of $`p_{MF}(a)`$. What the fluctuations do is that they provide a long range algebraic tail to the interactions, which overwhelms the exponentially screened mean-field attraction at large intersurface separation. The algebraic fluctuation tail decays with an anomalous exponent of $`1.7`$, that makes the polyelectrolyte mediated fluctuation interactions even longer ranged than the standard van der Waals interactions. This in itself is one of the main conclusions of our work. It is also this long ranged tail that is probably most amenable to experimental observation and verification. Since the fluctuation effects are stronger for smaller $`\alpha _{MF}`$ we predict that they should be easier to measure for very flexible polyelectrolytes, such as hyaluronic acid as opposed to DNA.

The main drawback of this work is that the fluctuation effects have been treated on a harmonic level which a priori assumes that they are small. However what saves us here is that the mean-field interactions decay exponentially whereas the fluctuation interactions decay algebraically. This is different than in the case of a simple, i.e. unconnected, Coulomb fluid such as an inhomogeneous electrolyte. In that case both the mean-field as well as the fluctuation interactions decay exponentially but the mean-field interactions have a range twice as long as the fluctuation interactions . The case of polyelectrolytes is in this respect fundamentally different: the range of fluctuation interactions is larger then the range of mean-field interactions and the use of harmonic approach to fluctuations is probably better grounded in this case.

## 7 Acknowledgements

The authors would like to thank R. Netz, P.L. Hansen and P. Ziherl for their comments on an earlier version of the manuscript. We also acknowledge the support of the Ministry of Science and Technology of Slovenia through the grant # J1-6168.


text:
# Localization and absence of Breit-Wigner form for Cauchy random band matrices

## Abstract

We analytically calculate the local density of states for Cauchy random band matrices with strongly fluctuating diagonal elements. The Breit-Wigner form for ordinary band matrices is replaced by a Levy distribution of index $`\mu =1/2`$ and the characteristic energy scale $`\alpha `$ is strongly enhanced as compared to the Breit-Wigner width. The unperturbed eigenstates decay according to the non-exponential law $`e^{\sqrt{\alpha t}}`$. We analytically determine the localization length by a new method to derive the supersymmetric non-linear $`\sigma `$ model for this type of band matrices.

Band matrices with random elements appear in a varity of physical problems in the context of quantum chaos and localization . A detailed analytical investigation of the properties of such matrices was performed by Fyodorov and Mirlin  using Efetov’s supersymmetry technique . Later, motivated by the localization problem of two interacting particles, Shepelyansky introduced and studied random band matrices superimposed with strongly fluctuating diagonal matrix elements . Subsequent work on this type of matrices  showed the existence of the Breit-Wigner regime where the eigenstates have a very peaked structure inside an eventual localization domaine.

In all these cases the matrix elements were drawn from a regular, typically gaussian, distribution with finite variance. In the present work, we extend the band matrix ensembles of refs.  to the case of Cauchy distributed matrix elements. In an early work, Lloyd  already introduced and studied a model with diagonal Cauchy distributed disorder. A detailed study of the level statistics of full random matrices  with elements drawn from the more general Levy distribution  was performed by Cizeau and Bouchaud. Band matrix ensembles with Cauchy distributed matrix elements were recently argued to be relevant and studied numerically in the context of the localization problem of two interacting particles . In this work, we present analytical results concerning the wave function properties for a model similar as in ref.  but which differs from that of ref.  by two important features: it concerns a banded and not full matrix, and its diagonal elements are typically much larger than the off-diagonal coupling elements.

Actually, we think that this case is of particular interest for a generic disordered fermionic many particle problem (for instance recently studied in ) where non-interacting eigenstates are coupled by quasi-random two-particle interaction matrix elements. There is an interesting and subtle regime  where the effective two-particle level spacing is smaller than the typical interaction matrix element while typical eigenstates are nevertheless composed of many non-interacting eigenstates. In this regime, we can diagonalize the problem for two particles by perturbation theory and apply the corresponding unitary transformation to the full many-particle Hamiltonian. This provides a new type of random matrix with residual three-body interaction matrix elements typically given by perturbative expressions with denominators, containing independent random non-interacting energy differences. Iterating this procedure one also obtains higher order contributions. Flambaum et al.  have argued that such denominators give effective distributions for these matrix elements characterized by long tails with the same power law as the Cauchy distribution. Even though the real situation is more complicated, this argument indeed shows the relevance of Cauchy random matrix ensembles.

The model we consider is a random band matrix of dimension $`N\times N`$ and of width $`b1`$ whose matrix elements are of the form $`H_{kl}=\eta _k\delta _{kl}+U_{kl}`$ with $`U_{kl}=0`$ if $`k=l`$ or $`|kl|>b`$. These matrix elements are statistically independent and distributed according to a Cauchy distribution $`p_a(x)\pi ^1a/(a^2+x^2)`$ where the width $`a`$ is given in terms of two parameters $`W`$ and $`U_0`$ via $`aW/\pi `$ for $`\eta _k`$ and $`aU_0`$ for $`U_{kl}`$. In this paper, we consider the orthogonal symmetry class where $`H`$ is real symmetric. The generalization to the other symmetry classes is straightforward. For the following analytical calculations, we will furthermore concentrate on the case where $`U_0`$ is sufficiently large to avoid a simple perturbative situation but still so small that the total density of states will essentially be dominated by the unperturbed diagonal elements only.

For ordinary band matrices  with a finite variance of $`U_{jk}`$, this case corresponds to the Breit-Wigner regime where the eigenstates are in addition to an eventual space localization also localized in the (unperturbed) energy space, i. e. only the sites $`j`$ with $`|E\eta _j|\mathrm{\Gamma }`$ essentially contribute to an eigenstate of energy $`E`$ where $`\mathrm{\Gamma }2\pi (2b)|U_{jk}|^2/WW`$ is the Breit-Wigner width. This can be seen in the local density of states at site $`j`$ which is a Lorentzian in $`(E\eta _j)`$ of width $`\mathrm{\Gamma }/2`$ .

To study the Cauchy band matrix ensemble introduced above, we therefore start with this quantity,

$$\rho _j(E)=\frac{1}{\pi }\text{Im}G_{jj}^{(+)}(E)_j,$$
(1)
where $`G^{(\pm )}(E)=(E\pm i0H)^1`$ and the average $`\mathrm{}_j`$ is taken with respect to all random matrix elements except the diagonal element $`\eta _j`$ at the site $`j`$ under consideration. As in the original Lloyd model  the average over the other diagonal elements $`\eta _k`$, $`kj`$ can be exactly performed by replacing in the definition of the Green function these elements according to: $`\eta _k+i\mathrm{\hspace{0.17em}0}iW/\pi `$. Using a general algebraic identity for the block inverse of a matrix, we obtain

$`\rho _j(E)`$ $`=`$ $`{\displaystyle \frac{1}{\pi }}\text{Im}{\displaystyle \frac{1}{E+i0\eta _j+i\mathrm{\Gamma }_j(U)/2}}_U,`$ (2)
$`\mathrm{\Gamma }_j(U)`$ $`=`$ $`2i{\displaystyle \underset{k,l(j)}{}}U_{jk}\stackrel{~}{G}_{kl}^{(+)}U_{lj},`$ (3)
where $`\mathrm{}_U`$ denotes the average over $`U_{jk}`$ only, and $`\stackrel{~}{G}^{(+)}=(E+iW/\pi \stackrel{~}{U})^1`$ with $`\stackrel{~}{U}`$ being the $`(N1)\times (N1)`$-matrix obtained from $`U`$ after elimination of the $`j`$-th row and $`j`$-th column. Eq. (2) reproduces the Lorentzian Breit-Wigner form of width Re$`[\mathrm{\Gamma }_j(U)]/2`$ for a given realization of the coupling matrix $`U`$. For ordinary band matrices this width is a self-averaging quantity and the further $`U`$-average does not modify the Breit-Wigner form. However, for the case of Cauchy band matrices $`\mathrm{\Gamma }_j(U)`$ is strongly fluctuating and the $`U`$-average will considerably modify the Breit-Wigner form. Now, we consider the limit $`|E|W`$ and we neglect $`\stackrel{~}{U}`$ in the definition of $`\stackrel{~}{G}^{(+)}`$ which is possible for the regime we consider. Using, $`\stackrel{~}{G}^{(+)}(i\pi /W)𝟙_{}`$, we obtain

$$\rho _j(E)=\frac{1}{\pi }\text{Im}\left\{i_0^{\mathrm{}}𝑑te^{it(E\eta _j)}\left[\varphi \left(\frac{2\pi U_0^2t}{W}\right)\right]^{2b}\right\},$$
(4)
with the function $`\varphi (z)`$ defined by

$$\varphi (z)\frac{1}{\pi }_{\mathrm{}}^{\mathrm{}}𝑑u\frac{1}{1+u^2}e^{zu^2/2},z0.$$
(5)
In the limit $`b1`$, the integral (4) is dominated by the behavior at small $`t`$. Using the approximation $`\mathrm{ln}\varphi (z)\sqrt{\frac{2}{\pi }z}+𝒪(z)`$, we obtain

$`\rho _j(E)`$ $`=`$ $`{\displaystyle \frac{1}{\alpha }}L_{1/2}\left({\displaystyle \frac{E\eta _j}{\alpha }}\right),\alpha ={\displaystyle \frac{16U_0^2b^2}{W}},`$ (6)
$`L_\mu (s)`$ $`=`$ $`{\displaystyle \frac{1}{2\pi }}{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑te^{its|t|^\mu }.`$ (7)
Here $`L_\mu (s)`$ represents a Levy distribution of index $`\mu `$  with the behavior $`L_\mu (s)s^{1\mu }`$ for $`|s|1`$. The expressions (6) and (7) provide the first important results of this work. They show that the local density of states is still a peaked function of $`(E\eta _j)`$ but with two important modifications. First, the Lorentzian Breit-Wigner form is replaced by the Levy distribution $`L_{1/2}`$, and second, the characteristic energy scale behaves as $`\alpha (4/\pi )b\mathrm{\Gamma }\mathrm{\Gamma }`$ where $`\mathrm{\Gamma }`$ is the Breit-Wigner width for ordinary band matrices (with a finite variance $`U_{jk}^2U_0^2`$). We mention that the approximation to neglect the contribution of $`\stackrel{~}{U}`$ in $`\stackrel{~}{G}^{(+)}`$ is valid for $`\alpha W`$, i. e. $`U_0bW`$.

The first modification has an important implication for the time-evolution of a state $`|\psi (t)>`$ initially localized at one site, $`|\psi (0)>=|j>`$. Then, the average of the amplitude $`<\psi (t)|j>`$ obeys a non-exponential decay law  of the form $`e^{\sqrt{\alpha t}}`$ instead of $`e^{\mathrm{\Gamma }t/2}`$. This type of behavior was for instance recently found in the context of many-body effects in cold Rydberg atoms . We conclude that it is a very general feature due to the long tail distribution of residual interaction matrix elements in many-body problems .

The second modification concerning the large enhancement of the characteristic energy scale can be qualitatively understood by the fact that for a given realization of $`U`$ the sum in eq. (3) is essentially determined by the typical maximal value $`2bU_0`$ of the $`2b`$ different matrix elements $`U_{jk}`$.

It is pretty obvious to generalize eqs. (6) and (7) to the case where the matrix elements $`U_{jk}`$ are drawn from the more general Levy distribution $`U_0^1L_\mu (U_{jk}/U_0)`$. In this case, we have to replace in (6) $`L_{1/2}`$ by $`L_{\mu /2}`$ and the characteristic energy scale is enhanced according to: $`\alpha b^{2/\mu 1}\mathrm{\Gamma }`$.

We now turn to a field theoretical formulation of the Cauchy band matrix model in terms of a supersymmetric non-linear $`\sigma `$ model . This formulation provides the access to the correlators of an advanced Green function $`G^{(+)}`$ and a retarded Green function $`G^{()}`$ which are for example useful to study the transport properties, the level statistics or the appearance of dynamical localization. For this, typically, the ensemble average of a generating functional of the form

$$F(J)=\text{Sdet}^{1/2}\left(E+(\frac{\omega }{2}+i\epsilon )\mathrm{\Lambda }H𝟙_\mathrm{𝟠}+𝕁\right)$$
(8)
is considered. For the details concerning the supersymmetry method and the notations, we refer to the standard literature . Here, we only remind that $`\mathrm{\Lambda }=𝟙_{}\text{diag}(\mathrm{𝟙},\mathrm{𝟙},\mathrm{𝟙},\mathrm{𝟙},\mathrm{𝟙},\mathrm{𝟙},\mathrm{𝟙},\mathrm{𝟙})`$ is the $`8N\times 8N`$ supermatrix whose eigenvalues distinguish between the matrix blocks associated to the advanced and retarded Green functions and that $`J`$ is a source matrix which is used to generate the Green function correlators by taking suitable derivatives with respect to its matrix elements .

Normally, $`H`$ contains a certain number of gaussian random variables and it is possible to perform the ensemble average analytically which provides after a some transformations and approximations the non-linear supermatrix model . For the case of Cauchy band matrices, this approach has to be modified because of the non-gaussian distribution of its matrix elements. Our strategy will be to perform analytically by a new method the ensemble average with respect to the Cauchy distributed diagonal matrix elements $`\eta _j`$ while the coupling matrix $`U`$ is kept fixed. This will provide a $`\sigma `$ model formulation for an arbitrary coupling matrix $`U`$. The ensemble average with respect to $`U`$ can then be performed at a later stage when possible and convenient. As a first application of this procedure, we will derive an analytic expression for the localization length of the Cauchy band matrix model.

The main idea is to replace the diagonal matrix elements $`\eta _j`$ by random phases $`e^{i2\phi _j}`$ via $`\eta _j=(W/\pi )\mathrm{tan}\phi _j`$. These phases are uniformly distributed on the unit cercle and can be formally identified with the scattering matrix of a chaotic cavity ideally coupled to one scattering channel. Describing the Hamiltonian of the chaotic cavity by an $`M\times M`$ random matrix $`h_j`$ drawn from the gaussian orthogonal ensemble (GOE) one can apply the following substitution  :

$$\eta _j=(W/\pi )\mathrm{tan}\phi _j=(W/\pi )A_j^{}h_j^1A_j$$
(9)
where $`A_j`$ is an $`M`$-dimensional vector with components $`(A_j)_l=\delta _{1,l}`$ and $`h_j`$ is a real symmetric gaussian random matrix with variance $`(h_j)_{lk}^2=(1+\delta _{lk})/M`$. The substitution (9) is valid in the limit $`M\mathrm{}`$ which will be taken at the end of the calculations. We mention that this approach is related to recent work concerning the derivation of the $`\sigma `$ model by performing phase averages in models involving unitary quantum maps for classically chaotic systems such as the kicked rotator  or rough billiards . The replacement (9) and the following steps follow very closely the method described in ref. .

Inserting (9) in the generating functional (8) and applying suitable transformations inside the superdeterminant, one obtains an effective model defined on a Hilbert space of dimension $`NM`$ and characterized by an effective Hamiltonian containing the GOE-matrices $`h_j`$ in its diagonal blocks and $`(W/\pi )A_j(D^1)_{jk}A_k^{}`$ in its off-diagonal $`jk`$-block. Here $`D`$ is the supermatrix appearing inside the superdeterminant of (8) with $`\eta _j`$ put to zero. This model is quite similar to the IWZ-model used in ref.  and following the standard procedure to derive the $`\sigma `$ model, we obtain $`F(J)_\eta =DQe^{(Q)}`$ with the action

$`(Q)`$ $`=`$ $`{\displaystyle \frac{1}{2}}\text{Str}\mathrm{ln}\left(E+(\frac{\omega }{2}+i\epsilon )\mathrm{\Lambda }U+i\frac{W}{\pi }Q+J\right),`$ (10)
where $`Q`$ contains in its diagonal blocks $`8\times 8`$ supermatrices $`Q_j`$ belonging to Efetov’s coset space for the orthogonal symmetry class . We mention some important points concerning this result : (i) it is exact since the saddle-point approximation to derive the $`\sigma `$ model becomes exact in the limit $`M\mathrm{}`$; (ii) its range of applicability goes far beyond the particular model considered here since it is valid for an arbitrary coupling matrix $`U`$ and its derivation only needs the Cauchy distribution for the diagonal matrix elements; (iii) it generalizes the above mentioned replacement rule for one-point functions to the case of two-point functions : $`\eta _k+i0\mathrm{\Lambda }i(W/\pi )Q_k`$ where $`Q_k`$ is now a dynamical variable over which we have to integrate.

To analyze the localization properties, we may put $`\omega =0`$ and expand the action (10) for the long wave length limit. In the regime $`U_0bW`$, this gives apart from the source term the action

$$(Q)\frac{1}{16}𝑑x\xi _x(U)\text{Str}[(_xQ(x))^2]$$
(11)
with the quantity $`\xi _x(U)(4\pi ^2/W^2)_{j=1}^bj^2U_{x,x+j}^2`$. If $`\xi _x(U)`$ were independent of $`x`$ the action (11) would correspond to the standard $`\sigma `$ model for quasi one-dimensional wires with the localization length $`\xi \xi _x(U)`$ . This is actually the case for gaussian distributed $`U_{xy}`$ for which $`\xi _x(U)\xi _x(U)_U=4\pi ^2b^3U_0^2/3W^2`$ is a selfaveraging quantity. However, for the Cauchy band matrix model $`\xi _x(U)`$ is strongly fluctuating and its average does not even exist. It is therefore necessary to take the $`x`$-dependence properly into account. A similar situation has recently been encountered by Rupp et al. who derived the 1d $`\sigma `$ model with an $`x`$-dependent “diffusion constant” from a hierachical random-matrix model for many-body states . According to this work all quantities that can be extracted from the 1d $`\sigma `$ model only depend on the rescaled length variable $`s(x)=_0^x𝑑y\xi _y^1(U)`$. This justifies the intuitive average procedure $`\xi ^1lim_L\mathrm{}\frac{1}{L}_0^L𝑑x\xi _x^1(U)\xi _x^1(U)`$. Fortunately, this average is finite and well defined for the Cauchy band matrix case (if $`b3`$). It can be easily evaluated by the same techniques used above for the local density of states \[see eqs. (2)-(7)\]. For $`b1`$, we obtain :

$$\xi =\frac{2\pi b^4U_0^2}{W^2}=b\gamma ,\gamma \frac{\pi }{16}\left(\frac{\alpha }{\mathrm{\Delta }}\right)$$
(12)
where the dimensionless parameter $`\gamma `$ counts the number of well coupled levels inside a strip of typical length $`b`$, $`\mathrm{\Delta }=W/2b`$ is the effective level spacing of such a strip, and $`\alpha `$ is the characteristic energy scale for the local density of states (6). (For the more general case of Levy distributed $`U_{jk}`$, we find: $`\xi b^{2+2/\mu }U_0^2/W^2`$.) The relation (12) compares to $`\xi b(\mathrm{\Gamma }/\mathrm{\Delta })`$ for the case of ordinary band matrices . However, it contradicts the expression $`\xi b\sqrt{\gamma }`$ which was with proper translation of notations numerically obtained from a similar Cauchy band matrix model in ref. . We attribute this to the extreme numerical difficulty to access the regime $`1\gamma b`$ where we expect (12) to be valid. For $`\gamma 1`$ we enter the perturbative regime with $`\xi b/\mathrm{ln}(\gamma ^1)`$ while for $`\gamma b`$ the localization length is likely to saturate at $`\xi b^2`$.

To verify this picture, we have therefore numerically determined the localization length by the recursive Green function method for $`0.001\gamma 300`$ and $`10b100`$. The inset of Fig. 1 shows the ratio $`\xi /b`$ as a function of $`\gamma `$. At first sight, it is indeed very difficult to verify the analytical result (12) and a direct numerical fit (e. g. for $`b=100`$ and $`2\gamma 30`$) gives the power law $`\xi /b\gamma ^{0.6}`$ rather than $`\xi /b\gamma `$. To understand this, we note that the number of well mixed levels inside one bandwidth obviously saturates at $`b`$ for $`\gamma >b`$. It appears therefore reasonable to replace in (12) $`\gamma `$ by an effective value $`\gamma _{\mathrm{eff}.}`$ that interpolates between $`\gamma `$ and $`b`$. Actually, we find that our numerical data can be well approximated by the scaling curve (see Fig. 1)

$$\frac{\xi }{b}\frac{1.5}{\mathrm{ln}(1+1.5/\gamma _{\mathrm{eff}.})},\gamma _{\mathrm{eff}.}\frac{\gamma }{[1+\sqrt{3.2(\gamma /b)}]^2}$$
(13)
which reproduces the perturbative limit for $`\gamma 1`$, the analytical result (12) for $`1\gamma b`$ and the behavior $`\xi b^2`$ for $`\gamma b`$. Note that eq. (13) provides for the case $`\gamma /b1`$ relative large corrections of order $`\sqrt{\gamma /b}`$ explaining the numerical difficulties to clearly identify this regime.

In summary, we have obtained analytical results for the local density of states and the localization length for Cauchy random band matrices. We find that the Breit-Wigner form is replaced by a Levy distribution with index $`\mu =1/2`$ and that the Breit-Wigner width $`\mathrm{\Gamma }`$ for ordinary band matrices is replaced by a new enhanced energy scale $`\alpha b\mathrm{\Gamma }`$. This energy scale also determines the localization length $`\xi b(\alpha /\mathrm{\Delta })`$. From the technical point of view, we have derived a $`\sigma `$ model formulation for arbitrary matrix ensembles with Cauchy distributed diagonal disorder.

The author thanks G. Caldara, B. Georgeot, J. Lages and D. Shepelyansky for useful discussions.