text:
# I. Introduction. Adler-Moser Polynomials

## I. Introduction. Adler-Moser Polynomials

In their 1929 paper J.L.Burchnall and T.W. Chaundy  examined the following (apparently elementary) question: What condition must be satisfied by two polynomials $`p(z)`$ and $`q(z)`$ in one variable $`z`$ in order the indefinite integrals

$$\left(\frac{p(z)}{q(z)}\right)^2𝑑z,\left(\frac{q(z)}{p(z)}\right)^2𝑑z$$
(1)
may be rational, i.e. expressible without logarithms. Under the assumption that $`p`$ and $`q`$ have no common or repeated factors the problem is reduced to finding polynomial solutions of the following differential-recurrence equation

$$p^{\prime \prime }q2p^{}q^{}+pq^{\prime \prime }=0$$
(2)
It was shown that solutions to (2) are Adler-Moser polynomials :

$$p=\theta _i,q=\theta _{i+1},\theta _0(z)=1,\theta _1(z)=z,\mathrm{deg}\theta _i=\frac{i(i+1)}{2}$$
(3)
The recurrence relation (2) for the Adler-Moser polynomials (8) can be rewritten as the first-order (three term) relation ,

$$\theta _{n+1}^{}\theta _{n1}\theta _{n+1}\theta _{n1}^{}=(2n+1)\theta _n^2$$
(4)
Considering it as a first order linear differential equation for $`\theta _{n+1}`$ we obtain recursively all solutions introducing an integration constant at each step. Thus the $`n`$th polynomial depends on $`n`$ free parameters

$$\theta _n=\theta _n(z;t_1,t_2,\mathrm{},t_n)$$

Polynomiality of all solutions is assured by the rationality of (1).

Several first examples of the Adler-Moser polynomials are:

$$\begin{array}{c}\theta _0=1\hfill \\ \theta _1=z\hfill \\ \theta _2=z^3+t_2\hfill \\ \theta _3=z^6+5t_2z^3+t_3z5t_2^2\hfill \\ \mathrm{}\hfill \end{array}$$
(5)
In (5), we set $`t_1=0`$ , since $`t_1`$ can be absorbed by the translation $`zz+t_1`$

It was observed in  that the roots of the consecutive polynomials $`\theta _i,\theta _{i+1}`$ are the equilibrium coordinates of $`i(i+1)/2`$ positive and $`(i+1)(i+2)/2`$ negative Coulomb charges (with values $`\pm 1`$ correspondingly) on the plane. Namely, the function

$$E=\underset{i<j=1}{\overset{n(n+1)/2}{}}\mathrm{ln}|x_ix_j|+\underset{i<j=1}{\overset{(n+1)(n+2)/2}{}}\mathrm{ln}|y_iy_j|\underset{i=1}{\overset{n(n+1)/2}{}}\underset{j=1}{\overset{(n+1)(n+2)/2}{}}\mathrm{ln}|x_iy_j|$$

has a critical point when $`x_i,y_i`$ are roots of $`\theta _n,\theta _{n+1}`$ respectively.

The constructive way to find the Adler-Moser polynomials is to use their determinant representation built by means of the Darboux-Crum transformations of the operator $`d^2/dz^2`$.

This construction can be demonstrated in the following diagram:

$$H_i=L_iM_i=M_{i+1}L_{i+1}H_{i+1}=L_{i+1}M_{i+1}=M_{i+2}L_{i+2}H_{i+2}\mathrm{}$$
(6)
where

$$L_i=\frac{\theta _i(z)}{\theta _{i1}(z)}\frac{d}{dz}\frac{\theta _{i1}(z)}{\theta _i(z)},M_i=\frac{\theta _{i1}(z)}{\theta _i(z)}\frac{d}{dz}\frac{\theta _i(z)}{\theta _{i1}(z)}$$

and $`L_0=M_0=d/dz`$.

From (6), we derive the intertwining differential operators $`D_i,U_i`$

$$U_i\frac{d^2}{dz^2}=H_iU_i,\frac{d^2}{dz^2}D_i=D_iH_iD_i=M_i\mathrm{}M_1,U_i=L_1\mathrm{}L_i$$
(7)
connecting $`H_0=d^2/dz^2`$ with $`H_i`$.

This, together with the obvious property $`D_iU_i=\frac{d^{2i}}{dz^{2i}}`$ leads to the Wronskian representation for the Adler-Moser polynomials. In more details

$$\theta _n=\mathrm{const}W[\psi _1,\mathrm{},\psi _n],\psi _n^{\prime \prime }=\psi _{n1},\psi _1=z$$
(8)
where the $`n`$th polynomial is (up to a constant factor) the Wronskian $`W=detd^i\psi _j/dz^i_{0in1,1jn}`$.

The intertwining operators are also expressible in terms of Wronskians. For instance

$$U_i[\psi ]=W[\psi _1,\mathrm{},\psi _i,\psi ]/W[\psi _1,\mathrm{},\psi _i]$$
(9)
## II. Indefinite integrals related to equilibrium of charges

By analogy with the previous section we pose the following question: when does an equilibrium condition for charges of different sign is equivalent to rationality of some indefinite integrals?

Before addressing this question we write the equilibrium equation for two types of charges in terms of polynomials:

###### Lemma 1

Let

$$E=\underset{1i<jn+m}{}Q_iQ_j\mathrm{ln}|z_iz_j|$$
(10)
be a real function of $`n+m`$ complex variables $`z_i,i=1\mathrm{}n+m`$ and

$$Q_i,Q_i=\{\begin{array}{c}1,i=1\mathrm{}n\hfill \\ \mathrm{\Lambda },i=n+1\mathrm{}m+n\hfill \end{array}$$

Then $`E`$ has a critical point at

$$z_i=\{\begin{array}{c}x_i,i=1\mathrm{}n\hfill \\ y_{in},i=n+1\mathrm{}m+n\hfill \end{array}$$

iff

$$\{p,q\}_\mathrm{\Lambda }:=\frac{d^2p(z)}{dz^2}q(z)2\mathrm{\Lambda }\frac{dp(z)}{dz}\frac{dq(z)}{dz}+\mathrm{\Lambda }^2p(z)\frac{d^2q(z)}{dz^2}=0$$
(11)
where $`p,q`$ are the following polynomials

$$p(z)=\underset{i=1}{\overset{n}{}}(zx_i),q(z)=\underset{i=1}{\overset{m}{}}(zy_i)$$

Proof: Is a calculation. We use the partial-fraction decomposition of $`\frac{\{p,q\}_\mathrm{\Lambda }}{pq}`$ around singular points $`x_i`$ and $`y_i`$.

Now we are in position to find an equivalent condition for $`p`$ and $`q`$ in terms of rational integrals.

Indefinite integrals of the algebraic function do not contain logarithmic terms if the function does not have singularities of the type $`z^1`$. The following lemma shows when certain functions, related to solutions of (11) , are free of simple poles and their integrals are rational.

###### Lemma 2

Let $`p(z)=_{i=1}^n(zx_i),q(z)=_{i=1}^m(zy_i)`$ do not have multiple or common roots. Then residues of simple poles of $`q^{2\mathrm{\Lambda }}/p^2`$ and $`p^{2/\mathrm{\Lambda }}/q^2`$ vanish:

$$Res_{z=x_i}\frac{q(z)^{2\mathrm{\Lambda }}}{p(z)^2}=0,i=1\mathrm{}n,Res_{z=y_i}\frac{p(z)^{2/\mathrm{\Lambda }}}{q(z)^2}=0,i=1\mathrm{}m$$

iff $`p`$ and $`q`$ satisfy (11).

Both indefinite integrals

$$\frac{q(z)^{2\mathrm{\Lambda }}}{p(z)^2}𝑑z,\frac{p(z)^{2/\mathrm{\Lambda }}}{q(z)^2}𝑑z$$
(12)
are rational iff $`p,q`$ satisfy (11) with $`\mathrm{\Lambda }=1/2,1,2`$

Proof: Let us factor $`p(z)`$ as $`p(z)=(zx_i)P(z)`$. Then the condition $`Res_{z=x_i}\frac{q(z)^{2\mathrm{\Lambda }}}{p(z)^2}=0`$ implies that

$$\frac{d}{dz}\left(\frac{q(z)^{2\mathrm{\Lambda }}}{P(z)^2}\right)_{z=x_i}=0$$

Since $`P(x_i)=p^{}(x_i)`$, $`P^{}(x_i)=p^{\prime \prime }(x_i)/2`$ we get the following equation

$$p^{\prime \prime }(x_i)q(x_i)2\mathrm{\Lambda }p^{}(x_i)q^{}(x_i)=0$$

Thus, since $`x_i`$ are roots of $`p(z)`$

$$p^{\prime \prime }(z)q(z)2\mathrm{\Lambda }p^{}(z)q^{}(z)=a(z)p(z)$$
(13)
where $`a(z)`$ is a polynomial.

Similarly, we get from the condition $`Res_{z=y_i}\frac{p(z)^{2/\mathrm{\Lambda }}}{q(z)^2}=0`$ that

$$\mathrm{\Lambda }^2p(z)q^{\prime \prime }(z)2\mathrm{\Lambda }p^{}(z)q^{}(z)=b(z)q(z)$$
(14)
Adding $`\mathrm{\Lambda }^2p(z)q^{\prime \prime }(z)`$ to (13) and $`p^{\prime \prime }q(z)`$ to (14) we get

$$p^{\prime \prime }(z)q(z)2\mathrm{\Lambda }p^{}(z)q^{}(z)+\mathrm{\Lambda }^2p(z)q^{\prime \prime }(z)=A(z)p(z)=B(z)q(z)$$

where $`A(z)=a(z)+\mathrm{\Lambda }^2q^{\prime \prime }(z),B(z)=b(z)+p^{\prime \prime }(z)`$.

Since $`p(z)`$ and $`q(z)`$ do not have common roots it follows from

$$A(z)p(z)=B(z)q(z)$$

that $`A(z)=N(z)q(z),B(z)=N(z)p(z)`$, where $`N(z)`$ is a polynomial. Thus, we arrive at

$$p^{\prime \prime }(z)q(z)2\mathrm{\Lambda }p^{}(z)q^{}(z)+\mathrm{\Lambda }^2p(z)q^{\prime \prime }(z)=N(z)p(z)q(z)$$
(15)
The degree of the polynomial in the left hand side of (15) is at most $`\mathrm{deg}(q)+\mathrm{deg}(p)2`$ while the degree of the right hand side is at least $`\mathrm{deg}(q)+\mathrm{deg}(p)`$, unless $`N(z)=0`$. Therefore, $`N(z)=0`$ and we obtain (11). The choice $`\mathrm{\Lambda }=1/2,\mathrm{\Lambda }=1,\mathrm{\Lambda }=2`$ is essential for the rationality of (12) because the both exponents $`2\mathrm{\Lambda }`$ and $`2/\mathrm{\Lambda }`$ are integer. This completes the proof.

Since the cases $`\mathrm{\Lambda }=1/2`$ and $`\mathrm{\Lambda }=2`$ are equivalent (one is connected with the other by permutation of $`p`$ with $`q`$) we have two second order bilinear differential relations corresponding to two types of rational integrals

* Adler-Moser polynomials, $`\mathrm{\Lambda }=1`$ (see )

$$\begin{array}{c}\{p,q\}_1=\frac{d^2p(z)}{dz^2}q(z)2\frac{dp(z)}{dz}\frac{dq(z)}{dz}+p(z)\frac{d^2q(z)}{dz^2}=0\hfill \\ \\ \frac{q(z)^2}{p(z)^2}𝑑z,\frac{p(z)^2}{q(z)^2}𝑑z\hfill \end{array}$$
(16)
* $`\mathrm{\Lambda }=2`$,

$$\begin{array}{c}\{p,q\}_2=\frac{d^2p(z)}{dz^2}q(z)4\frac{dp(z)}{dz}\frac{dq(z)}{dz}+4p(z)\frac{d^2q(z)}{dz^2}=0\hfill \\ \\ \frac{q(z)^4}{p(z)^2}𝑑z,\frac{p(z)}{q(z)^2}𝑑z\hfill \end{array}$$
(17)

Equation $`\{p,q\}_2=0`$ has been mentioned in  in connection with the problem of equilibrium of the point vortices in two-dimensional hydrodynamics.

We devote the sequel of the work to investigation of the second case and its generalizations.

## III. $`\mathrm{\Lambda }=2`$

We begin by summing up what we have obtained so far.

Let $`p(z)`$ and $`q(z)`$ be polynomials of the complex variable $`z`$ which do not have common or multiple roots.

We study the following three problems:

* When the indefinite integrals

$$\frac{p}{q^2}𝑑z,\frac{q^4}{p^2}𝑑z$$
(18)
are rational ?
* When the system of $`n`$ and $`m`$ Coulomb charges of values $`1`$ and $`2`$ respectively can be in equilibrium in the complex plane and what are coordinates of the charges ?

In other words when does the energy function

$$E=\underset{i<j}{}Q_iQ_j\mathrm{log}|z_iz_j|,Q_i=\{\begin{array}{c}1,i=1..n\hfill \\ 2,i=n+1..n+m\hfill \end{array}$$
(19)
have a critical point and at which $`z_i,i=1..n+m`$ ?
* When does the equation

$$\{p,q\}_2=0$$
(20)
have polynomial solutions?

The answer to questions (18), (19) and (20) is given by

###### Proposition 1

Let $`i`$, and (modulo translation of $`z`$) $`p_1=z,p_0=q_0=1,q_1=z`$. Then all solutions of the recurrence relations

$$\begin{array}{c}q_{i+1}^{}q_iq_{i+1}q_i^{}=(3i+1)p_i\hfill \\ p_i^{}p_{i1}p_ip_{i1}^{}=(6i1)q_i^4\hfill \end{array}$$
(21)
are polynomials of the degrees

$$\mathrm{deg}(p_i)=i(3i+2),\mathrm{deg}(q_i)=i(3i1)/2,i$$

They provide all solutions of (20)

$$\{p_i,q_i\}_2=\{p_i,q_{i+1}\}_2=0,i$$

Integrals (18) are rational iff $`p`$ and $`q`$ are connected by the bilinear equation (20) i.e. iff

$$p=p_i,q=q_i,\mathrm{or}p=p_i,q=q_{i+1}$$

The energy function (19) has a critical point provided $`z_i`$ are zeros of above pairs of polynomials: zeros of $`p`$, and $`q`$ being positions of charges 1 and -2 respectively.

Proof:

* The equivalence between the bilinear equation (20) and rationality of integrals (18) is a corollary of Lemma 2.
* The equivalence between (20) and existence of critical points of the energy function (19) is a corollary of Lemma 1.
* Since $`p=z^n+\mathrm{},q=z^m+\mathrm{}`$, from the highest symbol of (20) we get the Diophantine equation connecting $`n`$ with $`m`$

$$(n2m)^2n+4m=0$$

Its solutions are

$$n=n_i,m=m_i,\mathrm{or}n=n_i,m=m_{i+1},i$$

where

$$n_i=i(3i+2),m_i=\frac{i(3i1)}{2}$$

Let $`\{p_i,q_i\}_2=0`$. Considering this equation as a second order differential equation with a solution $`q_i`$, by elementary methods we find that its second linearly independent solution is given by

$$q_{i+1}=(3i+1)q_i\frac{p_i}{q_i^2}𝑑z,\{p_i,q_{i+1}\}_2=0$$
(22)
and is polynomial by the rationality of (18). The degrees of polynomials are connected by

$$\mathrm{deg}(q_{i+1})+\mathrm{deg}(q_i)=\mathrm{deg}(p_i)+1$$

It is seen that this relation is satisfied if $`\mathrm{deg}(p_i)=n_i,\mathrm{deg}(q_i)=m_i`$.

The similar procedure holds if we fix $`q_i`$ and consider $`p_{i1}`$ and $`p_i`$ as linearly independent solutions of (20):

$$p_i=(6i1)p_{i1}\frac{q_i^4}{p_{i1}^2}𝑑z$$
(23)
By freedom in choosing linearly independent solutions of (20) we write analogs of (22) and (23) for decreasing indices

$$q_i=(3i+1)q_{i+1}\frac{p_i}{q_{i+1}^2}𝑑z,p_{i1}=(6i1)p_i\frac{q_i^4}{p_i^2}𝑑z$$
(24)
and we can generate $`p_i,q_i`$ by induction in either directions starting at some $`i`$.

Since $`n_i`$ and $`m_i`$ are strictly increasing for $`i0`$ (strictly decreasing for $`i0`$) and $`m_0=n_0=0`$, the induction terminates for these two branches at $`p_0=q_0=1`$.

Rewriting (22), (23) or (24) in the differential form, we get the first-order recurrence relations (21), which completes the proof.

Here we write down several examples of polynomials satisfying above conditions for the branch $`i0`$

$$\begin{array}{cc}q_0=1\hfill & p_0=1\hfill \\ q_1=z\hfill & p_1=z^5+t_1\hfill \\ q_2=z^5+\tau _2z4t_1\hfill & p_2=z^{16}+\mathrm{}\hfill \\ \mathrm{}\hfill & \mathrm{}\hfill \end{array}$$

and for the branch $`i0`$

$$\begin{array}{ccc}p_0=1\hfill & & q_0=1\hfill \\ p_1=z\hfill & & q_1=z^2+\tau _1\hfill \\ \begin{array}{c}p_2=z^8+\frac{28}{5}\tau _1z^6+14\tau _1^2z^4\hfill \\ +28\tau _1^3z^2+t_2z7\tau _1^4\hfill \end{array}\hfill & & \begin{array}{c}q_2=z^7+7\tau _1z^5+35\tau _1^2z^3+\tau _2z^235\tau _1^3z\hfill \\ +\tau _1\tau _2\frac{5}{2}t_2\hfill \end{array}\hfill \\ \mathrm{}\hfill & \mathrm{}\hfill & \end{array}$$

where $`t_i`$ and $`\tau _i`$ stand for arbitrary constants.

## IV. Intertwining and factorization

One can try to find polynomials for $`\mathrm{\Lambda }=2`$ explicitly by analogy with the Adler-Moser case.

It easy to observe that equation (20) can be written in two different Scrödinger forms

$$\begin{array}{cc}S\left[\frac{p}{q^2}\right]=\left(\frac{d^2}{dz^2}+6(\mathrm{ln}q)^{\prime \prime }\right)\left[\frac{p}{q^2}\right]=0\hfill & (\mathrm{a})\hfill \\ \stackrel{~}{S}\left[\frac{q}{\sqrt{p}}\right]=\left(\frac{d^2}{dz^2}+\frac{3}{4}(\mathrm{ln}p)^{\prime \prime }\right)\left[\frac{q}{\sqrt{p}}\right]=0\hfill & (\mathrm{b})\hfill \end{array}$$
(25)
connected by the following factorization

$$\frac{p}{q^2}S=LM,\frac{p}{q^2}\stackrel{~}{S}=ML,L=\frac{d}{dz}\frac{\sqrt{p}}{q},M=\frac{p^{3/2}}{q^3}\frac{d}{dz}\frac{q^2}{p}$$

Since $`p=p_{i1}`$ and $`p=p_i`$ are solutions of (20) with $`q=q_i`$, we can put any of them in (25 a) with $`q=q_i`$. The similar fact holds for (25 b), with $`q=q_i`$ or $`q=q_{i+1}`$ and $`p=p_i`$. This ambiguity in the choice of solutions results in a possibility of different factorizations for the same second order operator. We remind that a similar possibility led to explicit representation of the Adler-Moser polynomials for $`\mathrm{\Lambda }=1`$.

Now, the factorization diagram has the following form

$$\begin{array}{cccccc}& & & & & \\ \frac{p_{i1}}{q_i^2}S_i=\stackrel{~}{L}_i\stackrel{~}{M}_i& & \frac{p_i}{q_i^2}S_i=L_iM_i& & & \\ & & & & & \\ & & \frac{p_i}{q_i^2}\stackrel{~}{S}_i=M_iL_i& & \frac{p_i}{q_{i+1}^2}\stackrel{~}{S}_i=\stackrel{~}{M}_{i+1}\stackrel{~}{L}_{i+1}& \\ & & & & & \\ & & & & \frac{p_i}{q_{i+1}^2}S_{i+1}=\stackrel{~}{L}_{i+1}\stackrel{~}{M}_{i+1}& \end{array}$$
(26)
where

$$\stackrel{~}{L}_i=\frac{d}{dz}\frac{\sqrt{p_{i1}}}{q_i},\stackrel{~}{M}_i=\frac{p_{i1}^{3/2}}{q_i^3}\frac{d}{dz}\frac{q_i^2}{p_{i1}},L_i=\frac{d}{dz}\frac{\sqrt{p_i}}{q_i},M_i=\frac{p_i^{3/2}}{q_i^3}\frac{d}{dz}\frac{q_i^2}{p_i}$$

The horizontal arrows in (26) correspond to factorizations connected with different choices of $`p`$ at fixed $`q`$ and vice versa, while the vertical arrows correspond to permutation of factors. Due to the change of (prefactors of) the second order operators in the horizontal lines we do not have simple intertwining relations between $`S_0=d^2/dz^2`$ and $`S_i`$:

$$U_i^LS_0=S_iU_i^R$$

Note that now there are two different $`U^L,U^R`$ (left, right) intertwining operators instead of one. This makes the line of approach used for Adler-Moser polynomials unpromising for $`\mathrm{\Lambda }=2`$.

## V. Baker-Akhieser functions, charges in homogeneous field

Again we ask the question: When the following indefinite integrals

$$\mathrm{exp}(kz)\frac{p(k,z)}{q(k,z)^2}\mathrm{exp}(kz)𝑑z,\mathrm{exp}(2kz)\frac{q(k,z)^4}{p(k,z)^2}\mathrm{exp}(2kz)𝑑z$$
(27)
are rational? In (27), $`k`$ is any complex number.

We call the function

$$\mathrm{\Psi }(k,z)=\frac{p(k,z)}{q(k,z)^2}\mathrm{exp}(kz)$$
(28)
rational Backer-Akhieser function for $`\mathrm{\Lambda }=2`$ by analogy with $`\mathrm{\Lambda }=1`$ (see below).

Similarly to the section III we can formulate the following

###### Proposition 2

Let polynomials $`p(k,z),q(k,z)`$ in $`z`$ do not have repeated or common roots. The following three statements are equivalent:

* The indefinite integrals

$$\mathrm{exp}(kz)\mathrm{\Psi }(k,z)𝑑z,\mathrm{exp}(2kz)\mathrm{\Psi }(k,z)^2𝑑z$$

are rational in $`z`$ functions, expressible without logarithms
* The energy function

$$\begin{array}{c}E=k_iQ_iz_i+_{i<j}Q_iQ_j\mathrm{ln}|z_iz_j|,\hfill \\ \\ Q_i=\{\begin{array}{c}1,i=1\mathrm{}n\hfill \\ 2,i=n+1\mathrm{}m+n\hfill \end{array},n=\mathrm{deg}(p),m=\mathrm{deg}(q)\hfill \end{array}$$
(29)
has a critical point at

$$z_i=\{\begin{array}{c}x_i(k),i=1\mathrm{}n\hfill \\ y_{in}(k),i=n+1\mathrm{}m+n\hfill \end{array}$$

where $`x_i,y_i`$ are roots of

$$p=\underset{i=1}{\overset{n}{}}\left(zx_i(k)\right),q=\underset{i=1}{\overset{m}{}}\left(zy_i(k)\right)$$
* $`p,q`$ satisfy the following equation

$$0=\{p,q\}_2^{(k)}:=\{p,q\}_2+2k(p^{}q2q^{}p)$$
(30)

Proof: Follows arguments similar to Lemmas 1, 2.

One can easily see that equation (29) describes equilibrium of two types of charges with values 1, -2 in the homogeneous electric field.

Another observation consists in the fact that

$$p(k,z)=p(\zeta ),q(k,z)=q(\zeta ),\zeta =kz$$
(31)
Indeed, one can easily check this fact by scaling $`z_i,i=1\mathrm{}n+m`$ in (29), or by change the integration variable $`z\zeta =kz`$ in (27).

Thus, we can set $`k=1`$ whenever $`k0`$.

Equating highest symbols in (30) we obtain the following simple

###### Lemma 3

Let $`p,q`$ satisfy conditions of proposition 2 then their degrees are related by

$$deg(p)=n=2m=2\mathrm{deg}(q)$$

In other words , the total charge has to be zero in order the system to be at rest in the homogeneous field.

Although the case $`\mathrm{\Lambda }=2`$ has some similarity with $`\mathrm{\Lambda }=1`$, the same approach to explicit representation of $`p_i,q_i`$ , as shown in the section IV, does not apply any more.

Here we write down several first examples of $`p_i(k,z),q_i(k,z)`$

$$q_0(k,z)=1,p_0(k,z)=1,$$

$$q_1(k,z)=\zeta ,p_1(k,z)=\zeta ^23\zeta +3,$$

$$\begin{array}{c}q_2(k,z)=\zeta ^3+t_2\zeta ^2+\frac{t_2^2+6}{3}\zeta ,\hfill \\ \\ p_2(k,z)=\zeta ^6+\left(9+2t_2\right)\zeta ^5+\left(4015t_2+5/3t_{2}^{}{}_{}{}^{2}\right)\zeta ^4\hfill \\ \\ +\left(96+52t_210t_{2}^{}{}_{}{}^{2}+2/3t_{2}^{}{}_{}{}^{3}\right)\zeta ^3+\left(11290t_2+\frac{76}{3}t_{2}^{}{}_{}{}^{2}3t_{2}^{}{}_{}{}^{3}+1/9t_{2}^{}{}_{}{}^{4}\right)\zeta ^2\hfill \\ \\ +\left(48+66t_228t_{2}^{}{}_{}{}^{2}+5t_{2}^{}{}_{}{}^{3}1/3t_{2}^{}{}_{}{}^{4}\right)\zeta 18t_2+10t_{2}^{}{}_{}{}^{2}3t_{2}^{}{}_{}{}^{3}+1/3t_{2}^{}{}_{}{}^{4}+48\hfill \end{array}$$
(32)
where $`t_2`$ is an arbitrary constant.

It is easy to show that solution (32) is generic for $`m=3`$ , in the sense that it is not contained in a wider class. In other words , coefficients of (properly normalized) $`q_2`$ cannot depend on more than one free parameter. Indeed , we can always take $`q_2`$ to be monic and without constant term (by fixing one root at $`z=0`$). If $`q_2`$ depended on more than one parameter, then it would be $`q_2=z^3+t_2z^2+Tz`$, where $`t_2`$ and $`T`$ are arbitrary constants. In such a case (30) become a linear equations for the coefficients of $`p_2`$. It is verified by elementary linear algebra that this system is incompatible. Thus $`q_2`$ can depend only on one free parameter (say $`t_2`$).

It is difficult to get examples of polynomials of higher degrees, even with the help of the computer. It is also nontrivial problem to define degrees of $`q(k,z)`$. If one assumes that by analogy with $`\mathrm{\Lambda }=1`$ case the order of $`q_i(k,z)`$ is a quadratic function of $`i`$ (which certainly holds for $`k=0`$, see proposition 1), then it follows from the above examples that $`\mathrm{deg}(q_i(k,z))=i(i+1)/2`$.

## VI. Scaling, reparametrization and bispectrality

A statement similar to proposition 2 holds for the case $`\mathrm{\Lambda }=1`$ with the Baker-Akhieser function $`\psi =(p/q)\mathrm{exp}(kz)`$ instead of (28) and symmetric integration conditions

$$\mathrm{exp}(kz)\psi 𝑑z,\mathrm{exp}(kz)\psi ^1𝑑z$$

The electrostatic analogy describes charges 1, -1 in the homogeneous field and equation (30) is replaced by

$$0=\{p,q\}_1+2k(p^{}qq^{}p)$$

This equation can be rewritten in the Schrödinger form

$$L\psi (k,z)=\left(\frac{d^2}{dz^2}+V\right)\psi (k,z)=k^2\psi (k,z),V=2\frac{d^2}{dz^2}\mathrm{log}q$$
(33)
It follows from (7), (9) (for more details see ) that the solution of (33) is obtained by action of the intertwining operator (9) to the eigenfunction $`\mathrm{exp}(kz)`$ of the free Schrodinger operator $`d^2/dz^2`$. This solution is the Backer-Akhieser function:

$$\psi =W[\psi _1,\mathrm{},\psi _n,\mathrm{exp}(kz)]/W[\psi _1,\mathrm{},\psi _n]$$

Therefore, according to (8)

$$p(k,z)=\mathrm{exp}(kz)W[\psi _1,\mathrm{},\psi _n,\mathrm{exp}(kz)],q(k,z)=W[\psi _1,\mathrm{},\psi _n]=\mathrm{const}\theta _n(z)$$

On the other hand , by arguments leading to (31) , we see that

$$\left(\frac{d^2}{d\zeta ^2}+V\right)\psi (\zeta )=\psi (\zeta ),V=2\frac{d^2}{d\zeta ^2}\mathrm{log}q$$
(34)
with

$$\psi (\zeta )=\frac{p(\zeta )}{q(\zeta )}\mathrm{exp}(\zeta ),p(\zeta )=\mathrm{exp}(\zeta )W[\psi _1(\zeta ),\mathrm{},\psi _n(\zeta ),\mathrm{exp}(\zeta )],q(\zeta )=\mathrm{const}\theta _n(\zeta )$$

We observe that potential $`V=V(z)`$ in (33) is independent of $`k`$, while potential $`V=V(\zeta )`$, in (34) , depends on $`k`$. This is due to the fact that rescaling of variable $`z\zeta =kz`$ can be absorbed by appropriate changes of the parameters (see e.g. (5))

$$\theta _n(z;t_1,t_2,\mathrm{},t_n)=k^{n(n+1)/2}\theta _n(kz;kt_1,k^3t_2,\mathrm{},k^{2n1}t_n)$$

since the Adler-Moser polynomials can be seen as homogeneous polynomials of variables $`z`$ and $`t_i,i>0`$ with weights 1 and $`2i1`$ correspondingly.

Thus the problem of equilibrium of charges $`1`$ and $`1`$ can be reduced to the spectral problem for the $`k`$-independent operator $`L`$ (33).

Similarly, one can ask if $`\psi (k,z)`$ is simultaneously a solution of a spectral problem for $`z`$ independent differential with respect to $`k`$ operator $`A`$

$$A\psi (k,z)=\mathrm{\Theta }(z)\psi (k,z),A/z=0$$

where $`\mathrm{\Theta }(z)`$ is a function of $`z`$.

It turns out  that such operators exist and belong to the so called odd bispectral family.

In general, the problem of finding functions satisfying simultaneously a differential equation in $`z`$ with $`k`$-dependent eigenvalues and a differential equation in $`k`$ with $`z`$-dependent eigenvalues is called the bispectral problem .

Now returning to the main subject of this work , one might expect that the $`\mathrm{\Lambda }=2`$ case can belong to some other bispectral (e.g. even ) family of the differential operators. Unfortunately , the $`\mathrm{\Lambda }=2`$ case is not , in general , related to the bispectral problem for a second order differential operator. Indeed , equation (30) can be rewritten in the following two forms

$$\begin{array}{c}\left(\frac{d^2}{d\zeta ^2}+u(\zeta )\right)\mathrm{\Psi }(\zeta )=\mathrm{\Psi }(\zeta ),u(\zeta )=6\frac{d^2}{d\zeta ^2}\mathrm{log}q(\zeta )\hfill \\ \\ \left(\frac{d^2}{d\zeta ^2}+v(\zeta )\right)\mathrm{\Psi }(\zeta )^{1/2}=\mathrm{\Psi }(\zeta )^{1/2},v(\zeta )=\frac{3}{4}\frac{d^2}{d\zeta ^2}\mathrm{log}p(\zeta )\hfill \end{array},\mathrm{\Psi }(\zeta )=\frac{p(\zeta )}{q(\zeta )^2}\mathrm{exp}(\zeta )$$
(35)
We use (32) as a counterexample to the bispectrality. It is seen that there is no change of parameter $`t_2f(t_1,t_2)`$ combined with the translation $`zz+t_1`$ in (32) , such that $`p_2`$ or $`q_2`$ become homogeneous (with some weights) in $`z`$ and $`t_2`$. Thus, in difference from the $`\mathrm{\Lambda }=1`$ case, rescaling of variable $`z`$ cannot be absorbed by a change of parameter $`t_2`$. As a consequence , potentials $`u,v`$ in (35) cannot be $`z`$ or $`k`$ independent. Therefore , it is impossible , in general , to set a bispectral problem with the second order differential operator in $`z`$ (or $`k`$) when $`\mathrm{\Lambda }=2`$.

## VII. Integrable dynamics of charges

The charge configurations studied above can be considered as fixed points of some dynamical systems .

###### Lemma 4

Let two polynomials $`p(z,t)=_{i=1}^n(zz_i(t))`$ and $`q(z,t)=_{i=n+1}^{n+m}(zz_i(t))`$ satisfy the following bilinear equation

$$q\frac{dp}{dt}\mathrm{\Lambda }p\frac{dq}{dt}=\{p,q\}_\mathrm{\Lambda }$$
(36)
Then the roots $`z_i(t),i=1\mathrm{}n+m`$ satisfy the system of ordinary differential equations

$$\frac{dz_i}{dt}=\underset{ji=1}{\overset{n+m}{}}\frac{Q_j}{z_iz_j}=\frac{1}{Q_i}\frac{E}{z_i}$$
(37)
where

$$Q_i=\{\begin{array}{c}1,i=1\mathrm{}n\hfill \\ \mathrm{\Lambda },i=n+1\mathrm{}m+n\hfill \end{array}$$

It follows from the above lemma that the critical points of the energy (10) are fixed points of (37). Although, as seen from the previous considerations, existence of such critical points depends on values $`m,n`$ and $`\mathrm{\Lambda }`$, equation (36) has polynomial solutions for any integer $`m,n0`$ and real $`\mathrm{\Lambda }`$.

Although (37) (or equivalently (36) ) is not a Hamiltonian system, it can be embedded in the Hamiltonian flow:

###### Proposition 3

(37) is a trajectory of a system with the Hamiltonian

$$H=\underset{i=1}{\overset{n+m}{}}\left(\frac{dz_i}{dt}\right)^2\underset{i<j=1}{\overset{n+m}{}}\frac{Q_iQ_j(Q_i+Q_j)}{(z_iz_j)^2}$$
(38)
In other words, the Hamiltonian equations of motion

$$\frac{d^2z_i}{dt^2}=\underset{ji=1}{\overset{n+m}{}}\frac{Q_iQ_j(Q_i+Q_j)}{(z_iz_j)^3}$$

are corollaries of (37)

Proof: This proposition is a corollary of more general result proved in  .

This is a nontrivial fact that the Hamiltonian system of the point masses interacting through two body forces is a corollary of equations of motion of a lower order for pairwise interacting points. In general , such a property is equivalent to compatibility of a highly overdetermined system of equations. Perhaps , elliptic generalizations of (37) are the only systems satisfying such conditions .

It is easy to note that, when $`\mathrm{\Lambda }=1`$ (38) is a sum of two independent Hamiltonians

$$H=H_1+H_2,H_1=\underset{i=1}{\overset{n}{}}\left(\frac{dz_i}{dt}\right)^2\underset{i<j=1}{\overset{n}{}}\frac{2}{(z_iz_j)^2},H_2=\underset{i=n+1}{\overset{n+m}{}}\left(\frac{dz_i}{dt}\right)^2\underset{i<j=n+1}{\overset{n+m}{}}\frac{2}{(z_iz_j)^2}$$

Each of them belongs to the completely integrable (in the Liouville sense) Calogero-Moser system . As a corollary (37) is integrable when $`\mathrm{\Lambda }=1`$ .

Equation (38) is also the Calogero Moser Hamiltonian if $`\mathrm{\Lambda }=1`$.

In general, (38) is not reduced to known integrable Hamiltonians.

It is conjectured in  that (37) is integrable for any real $`\mathrm{\Lambda }`$ in the sense that there exist $`2(n+m)1`$ functionally independent real constants of motion which are rational functions of $`z_i`$ and $`z_i^{}`$ (asterisk denotes the complex conjugation and $`i=1\mathrm{}n+m`$).

It can be seen from the above the case $`\mathrm{\Lambda }=1`$ has a particular significance both in the problem of the rational integrals (1) and from the point of view of the integrable Hamiltonian systems. It is interesting to establish the role of the second case $`\mathrm{\Lambda }=2`$ in the theory of the dynamical system (37).

## VIII Conclusion

The principal question of this paper to be addressed is about existence and explicit representation of the Baker-Akhieser function (or equivalently $`p(k,z),q(k,z)`$) for $`\mathrm{\Lambda }=2`$. It is worth to mention a useful generalization (interesting by itself) which might help to find the answer: By analogy with the $`\mathrm{\Lambda }=1`$ case, one may consider equilibrium of charges on the cylinder. Remind that (e.g. see ,), equilibrium configurations for $`\mathrm{\Lambda }=1`$ are zeros of trigonometric polynomials obtained by finite numbers of Darboux transformations from the free Schrödinger operator on circle. In such setting, both the Adler-Moser polynomials and $`\mathrm{\Lambda }=1`$ rational Baker-Akhieser function can be obtained choosing different parameters in the limit when the radius of the cylinder goes to infinity.

It might be interesting to examine a similar procedure getting an analog of proposition 1 on the cylinder.

Another interesting question , as was mentioned before , is to understand the role of the case $`\mathrm{\Lambda }=2`$ in the theory of dynamical systems (37) and integrability of related Hamiltonians.

Finally , we would like to mention relation with the Reduced Model of Superconductivity (Richardson, ) , the Gaudin Magnet Model  from the one hand and the $`\mathrm{\Lambda }=2`$ equilibrium configurations from the other hand. These models are sets of commuting quantum Hamiltonians acting on a finite-dimensional Hilbert space. Equations for the Hamiltonian eigenvalues are equilibrium conditions for a set of charges of value $`2`$ located in the complex plane, subject to mutual repulsion , and attraction of charges of value $`1`$ located at positions defined by the parameters of the Hamiltonian. The $`\mathrm{\Lambda }=2`$ case considered in the present work corresponds to the Richardson (Gaudin) models with special restrictions imposed on the parameters of Hamiltonian.It is interesting to examine relations between such models and establish their special properties connected with the $`\mathrm{\Lambda }=2`$ case.

## Acknowledgments

The author is grateful to H.Aref, Y.Berest, B.Dubrovin, T.Grava, O.Yermolayeva for useful information , help and remarks.


format: text

text:
# Epidemic processes with immunization

## I Introduction

Epidemic processes can be described as the spread and decay of a non-conserved agent, an example of which is an infectious disease Mollison77 . The agent is not allowed to appear sponstaneously but it can multiply itself by infecting neighbouring individuals, or decay at a constant rate. Depending on the balance between these two processes, the infection may either die out or spread over the entire population. The two regimes of survival and extinction of the epidemic are typically separated by a continuous non-equilibrium phase transition. When the decay process dominates, the epidemic dies out at large times and the system gets trapped in an absorbing state from which it cannot escape.

Continuous phase transitions into absorbing states are associated with certain universality classes MarroDickman ; Hinrichsen00 . For epidemic processes, a well studied case is the universality class of directed percolation (DP). It is believed that two-state spreading processes with short-range interactions generically belong to the DP class, provided that quenched randomness, unconventional symmetries and large scales due to memory effects are absent Janssen81 ; Grassberger82 . Examples of physical systems whose critical behaviour is described by DP include heterogeneous catalysis ZGB86 , chemical reactions GdelaT79 ; Schlogl72 , interface depinning TangLeschhorn92 ; Buldyrev , the onset of spatio-temporal chaos RRR03 , flowing sand FlowingSand  and self-organized criticality MohantyDhar .

The epidemic process in which the susceptibility to infection is independent of previous infections is described by DP. However, for a more realistic description, we should consider an immunization effect Perelson97 . Immunization can be added to the DP model by changing the susceptibility after the first infection Cardy83 ; CardyGrassberger85 ; Janssen85 . A minimalistic model that captures this feature is one that is controlled by two independent parameters: a probability of first infection and another probability for all subsequent reinfections. The fact that the local susceptibility depends on whether a site has been infected in the past or not leads to a non-Markovian epidemic process, in which the time evolution depends on the entire history. This non-Markovian feature changes the universality class of the epidemic spreading.

The phase diagram of an epidemic process with immunization (Fig. 1) was studied in Ref. GCR97 ; MGD98 . If the probabilities for first infections and reinfections are equal, the model corresponds to ordinary directed percolation. However if the susceptibility changes to zero after the first infection there is perfect immunization, and the model reduces to the General Epidemic Process (GEP) Mollison77 . GEP belongs to the ordinary percolation universality class Alexandrowicz80 . The critical points of the GEP and DP are connected by a curved phase transition line separating a phase in which the spreading process always dies out, from another phase of annular growth, where an active front may propagate into regions of non-immune sites, leaving a bulk of immune sites behind. As shown in Ref. GCR97  the critical behavior along this line (except for the upper terminal point) belongs to the same universality class as GEP. Using field theoretic renormalization group techniques, the critical exponents were calculated along this line CardyGrassberger85 ; Janssen85 .

The main result in Ref. GCR97  is that the compact growth/no growth phase transition line is at the critical value of the reinfection rate and independent of the first infection rate. Above this horizontal transition line in Fig. 1, the model exhibits compact growth and approaches the stationary state of supercritical DP. This is because, in the active phase of an epidemic process with immunization, each site will be visited at least once after a sufficiently long time so that the dynamics in the stationary active state involves only reinfections. On the horizontal phase transition line itself all reinfection processes are critical DP while the probability of first infections may be sub- or supercritical. By varying the first infection rate, we can impede or facilitate the spreading into non-immune regions. In Ref. GCR97  a numerical analysis of the scaling behaviour along this horizontal line gives no hint of power law behaviour, and it is also suggested that this result could be applied to models with multiple absorbing states MGD98 ; JensenDickman93 ; MDHM94 ; Dickman96 ; Munoz96 ; LopezMunoz97 ; Munoz97 ; HwangPark99 ; FvW02 .

The aim of the present work is to investigate in further detail the dynamical critical behavior along the horizontal phase transition line and in the vicinity of the DP critical point, and to give a theoretical explanation of the absence of power law scaling along this line.

In order to describe the effect of immunization in an epidemic process, in Sec. II we study a field theoretic formulation of directed percolation with immunization and we show that the non-Markovian term which contains the immunization effects, is relevant under renormalization group analysis. As a consequence of this result we argue that the asymptotic spreading properties along the horizontal transition line should be determined by the limits of very low and very high first infection probability. Therefore, in Sec. III we present a study of the very small infection probability limit and we develop a quasi-static approximation to obtain the scaling behaviour of the survival probability of the epidemic. It turns out that this does not follow a power law, but instead decays asymptotically as a stretched exponential. This theoretical prediction is confirmed by optimized numerical simulations in 1+1 dimensions. The high first infection probability limit is studied in Sec. IV giving similar results. We complete the section with a theoretical approximation for the spreading behaviour, and with numerical calculations to corroborate the theoretical claims. Finally, a discussion of these results together with a possible connection with multiple absorbing state models is the subject of the conclusions in Sec. V.

## II Field theoretical analysis

### II.1 The model

In this section we develop an alternative derivation of the action for the DP model with immunization that was proposed in Cardy83 ; CardyGrassberger85 ; Janssen85 . The microscopic rules for DP in $`d+1`$ dimensions are rather simple. An infected site at time $`t`$ can infect its nearest neighbours at time $`t+1`$ with a probability $`p_0`$. There is a critical threshold $`p_c`$ such that for $`p_0<p_c`$ the epidemic process always dies out, that is, it reaches the absorbing state. For $`p_0>p_c`$ there is a finite probability that the epidemics survives. At the critical point $`p_0=p_c`$ the system scales anisotropically in time and space. The upper critical dimension is $`d_c=4`$, below which the fluctuation effects become important. The field theoretic action of DP GrassbergerSundermeyer78 ; CardySugar80  reads as follows:

$$S_{\mathrm{DP}}=𝑑td^dx\left[\stackrel{~}{\varphi }(_tD^2+r)\varphi +u_1\stackrel{~}{\varphi }\varphi ^2u_2\stackrel{~}{\varphi }^2\varphi \right].$$
(1)
Here, $`\varphi `$ is the local activity, $`\stackrel{~}{\varphi }`$ is the response field and $`rp_cp_0`$ is the mass parameter which measures the distance from criticality. This action can also be written as a Langevin-type equation for the local activity,

$$(_tD^2+r)\varphi +\frac{1}{2}u\varphi ^2+\xi (𝐱,t)=0,$$
(2)
where $`u`$ is the symmetrized coupling constant after rescaling the fields according to the DP time reversal symmetry. The noise $`\xi (𝐱,t)`$ is Gaussian, and satisfies $`\xi (𝐱,t)=0`$ and $`\xi (𝐱,t)\xi (𝐱^{},t^{})=u\delta ^d(𝐱𝐱^{})\delta (tt^{})`$. For a systematic analysis of the immunization around the DP fixed point it is often more convenient to make use of the description in terms of the action.

To add immunization to the model, the susceptibility after the first infection is decreased by an amount $`\lambda >0`$. The probability of infection, $`p`$, depends locally on position and time, $`p=p(𝐱,t)`$. The microscopic rules are modified as follows: a healthy site can first be infected with infection probability $`p_0`$. A site which is infected at time $`t`$ will become immune at the next time step $`t+1`$. Meanwhile any immune site can be re-infected at the lower re-infection probability $`p_0\lambda `$. Thus, the state of the system at time $`t`$ depends not just on the configuration of infected sites at time $`t1`$, but on the entire previous history.

In order to take into account the effects of immunization, we will modify the DP action. Let us consider a discrete $`d+1`$ dimensional lattice and assume that the field $`\varphi (𝐱,t)`$ is almost constant between times $`t`$ and $`t+\mathrm{\Delta }t`$, where $`\mathrm{\Delta }t=1`$ is the time step unit on the lattice. We then subdivide it into $`N`$ intervals, $`\mathrm{\Delta }t/N`$. First of all, we want to find an expression for the probability that there is an active site in the temporal sub-interval $`\mathrm{\Delta }t/N`$. To do so in terms of the field $`\varphi (𝐱,t)`$, we need to take into account the fact that $`\varphi `$ has units of $`length^{d/2}`$. Then this probability can be written as $`w\varphi (𝐱,t)\mathrm{\Delta }t/N`$, where we introduce a parameter $`w`$ to ensure that this expression is dimensionless.

In this way, $`1w\varphi (𝐱,t)\frac{\mathrm{\Delta }t}{N}`$ is the probability that the site is not active in the interval $`\mathrm{\Delta }t/N`$. The probability of not finding activity between times $`t`$ and $`t+\mathrm{\Delta }t`$, can then be expressed as:

$`{\displaystyle \underset{i=0}{\overset{N}{}}}(1w\varphi (𝐱,t+i/N){\displaystyle \frac{\mathrm{\Delta }t}{N}})`$ $``$ $`(1w\varphi (𝐱,t){\displaystyle \frac{\mathrm{\Delta }t}{N}})^N`$ (3)
$`\stackrel{N\mathrm{}}{}`$ $`e^{w\varphi (𝐱,t)\mathrm{\Delta }t}.`$
However, between times $`t=0`$ and $`t>1`$, we can no longer assume that the field $`\varphi `$ is constant. Therefore, the probability of a site $`𝐱`$ has never been infected by time $`t`$ turns out to be,

$$\underset{t^{}<t}{\overset{t}{}}e^{w\varphi (𝐱,t^{})\mathrm{\Delta }t}=e^{w_{t^{}=0}^t\varphi (𝐱,t^{})\mathrm{\Delta }t}.$$
(4)
When the continuous limit is taken on the discrete lattice, the probability of having been infected at least once in the past becomes $`1\mathrm{exp}(w_0^t\varphi (𝐱,t^{})𝑑t^{})`$. This expression is the probability for a site to be immune at time $`t`$. Since $`\lambda `$ is the parameter related to the immunization in the system, we can use the above result to write an expression for $`p(𝐱,t)`$,

$$p(𝐱,t)=p_0\lambda \left[1e^{w_0^t\varphi (𝐱,t^{})𝑑t^{}}\right].$$
(5)
We assume that the mass parameter $`r`$ can be written as $`rp_cp(𝐱,t)`$. Thus, using Eq. 5, the addition of immunization can be reflected in the action as a modification in the original DP mass parameter by the following substitution:

$$rr+\lambda \left[1e^{w_0^t\varphi (𝐱,t^{})𝑑t^{}}\right].$$
(6)
Finally the action for DP with immunization can be written as a modification of the directed percolation action as follows:

$`S`$ $`=`$ $`S_{\mathrm{DP}}+`$ (7)
$`\lambda `$ $`{\displaystyle 𝑑td^dx\stackrel{~}{\varphi }(𝐱,t)\varphi (𝐱,t)\left[1e^{w_0^t𝑑t^{}\varphi (𝐱,t^{})}\right]}.`$
According to the field theory, $`w`$ is a finite coupling constant. Then in the numerical calculations on a discrete $`d+1`$ dimensional lattice, $`w`$ can be considered as a parameter such that the subsequent infection rate is,

$$p_0+(pp_0)\left[1\mathrm{exp}(wn(𝐱,t))\right].$$
(8)
Here, the function $`n`$ counts all the past activity at a site $`𝐱`$ until time $`t`$. Since the exponential function vanishes at sufficient large times, the effective subsequent infection rate is equal to $`p`$.

In the numerical calculations carried out in the present work, we assume that the susceptibility for spreading only changes at the first infection process and remains constant thereafter. This assumption is equivalent to take an infinite value of $`w`$, but does not make any relevant change in the final results of the theory as it is argued in Ref. GCR97 .

### II.2 Renormalization group analysis

Expanding in Eq. 7 the exponential function as a power series, the action reads

$`S`$ $`=`$ $`S_{\mathrm{DP}}+{\displaystyle }dtd^dx\stackrel{~}{\varphi }(𝐱,t)\varphi (𝐱,t)\times `$ (9)
$`\left[{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{\lambda ^{(n)}}{n!}}\left({\displaystyle _0^t}𝑑t^{}\varphi (𝐱,t^{})\right)^n\right],`$
where we abbreviate $`\lambda ^{(n)}=\lambda w^n`$. Let us now study the stability of the DP fixed point with respect to a small perturbation in $`\lambda ^{(n)}`$ by dimensional analysis. Introducing a length scale $`k^1`$, the mean-field fixed point of DP is characterized by the dimensions $`[\stackrel{~}{\varphi }\varphi ]=k^d`$, $`[D]=\omega k^2`$, and $`[u/D]=k^{4d/2}`$ so that the upper critical dimension is $`d_c=4`$. When immunization is considered as a perturbation, the effective expansion parameter in the expression of the two point vertex function is $`u_2^n\lambda ^{(n)}`$. Thus, by dimensional analysis we find

$$[\frac{u_2^n\lambda ^{(n)}}{D_0^{2n+1}}]=k^{2+n(4d)}.$$
(10)
Consequently the upper critical dimension for each coupling constant $`\lambda ^{(n)}`$ is $`d_c=\frac{2}{n}+4`$, where $`n=1,\mathrm{},\mathrm{}`$. This means that for dimensions $`d>6`$, $`\lambda ^{(n)}`$ is irrelevant for all $`n`$. The parameter $`u_2`$ will also be irrelevant, and only the mass parameter $`r`$ renormalizes.

Consider now the renormalization group (RG) parameter space spanned by $`r`$, $`u_2`$, and the $`\lambda ^{(n)}`$. The DP fixed point occurs at a point on the $`u_2`$ axis, and the GEP fixed point is somewhere in the hyper-dimensional space. Since all the coupling constants are irrelevant in the case $`d>6`$, the RG flows bring the system from the neighborhood of the DP fixed point towards the Gaussian DP fixed point, which corresponds to the DP free-field theory. Therefore if close to the DP fixed point we turn on a small perturbation of immunization in a system, the qualitative nature of the system will not change, and the epidemic process is expected to spread according to a mean-field DP process controlled by the probability of first infections.

For $`4<d6`$, the coupling constant $`\lambda ^{(1)}`$ is the most relevant. Consequently the higher order terms in the expansion of the exponential in Eq. (9) can be neglected. In this case the RG flows bring the system away from the neighbourhood of the DP fixed point to the GEP fixed point which is stable. Thus, we expect the system to undergo an ordinary crossover from DP to the dynamic percolation universality class at a certain typical time scale. However, recently it has been suggested that dangerous irrelevant operators may possibly lead to a non-trivial critical behavior different from dynamical percolation FvW02 .

To take the fluctuation effects into account, we define the renormalized $`\lambda ^{(1)}`$ as $`\lambda _R^{(1)}`$ and we apply standard methods of the perturbative renormalization group Amit78 . In particular we make use of an $`ϵ=4d`$ expansion around the DP critical point. The non-Markovian modification in the DP action can be written then as an additional term of the form,

$$\lambda ^{(1)}𝑑td^dx\mathrm{\Theta }(tt^{})\stackrel{~}{\varphi }(𝐱,t)\varphi (𝐱,t)\varphi (𝐱,t^{}).$$
(11)
The non-locality in time, expressed by the $`\mathrm{\Theta }`$ function in Eq. 11, is represented by a dashed line propagator in the Feynman diagrams. The renormalized coupling constant $`\lambda _R^{(1)}`$ can be determined in terms of the vertex function $`\mathrm{\Gamma }_{\stackrel{~}{\varphi }\varphi \varphi }^{}`$ (see Fig. 2), evaluated at some specified normalization point, NP, setting a momentum scale $`\mu `$. The  in the vertex function’s notation in Eq. 11 is used to indicate that this function is calculated in the DP with immunization theory and is different from the one calculated in the ordinary DP. Near the upper critical dimension the ultraviolet divergences are absorbed in the renormalization constant $`Z_\varphi `$, where $`\varphi _RZ_\varphi ^{1/2}\varphi `$ and $`Z_\varphi =Z_{\stackrel{~}{\varphi }}`$.

$`\lambda _R^{(1)}`$ $`=`$ $`\mathrm{\Gamma }_{R\stackrel{~}{\varphi }\varphi \varphi }^{}|_{NP}`$ (12)
$`=`$ $`Z_\varphi ^{3/2}\mathrm{\Gamma }_{\stackrel{~}{\varphi }\varphi \varphi }^{}|_{NP}.`$
The structure of the corrections to the vertex function $`\mathrm{\Gamma }_{\stackrel{~}{\varphi }\varphi \varphi }^{}`$ suggests a correspondence with the RG theory of DP away from the critical point, where is considered the renormalization of a nonzero mass term $`r`$ which couples with the composite operator $`\stackrel{~}{\varphi }(x,t)\varphi (x,t)`$. One loop corrections to the corresponding two point bare vertex function $`\mathrm{\Gamma }_{\stackrel{~}{\varphi }\varphi }`$ in this theory, are depicted in Fig. 3. The normalization condition implies that the renormalized two point vertex function $`\mathrm{\Gamma }_{R\stackrel{~}{\varphi }\varphi }`$ evaluated at the normalization point NP, is equal to 1.

$$\mathrm{\Gamma }_{R\stackrel{~}{\varphi }\varphi }|_{NP}=Z_{\stackrel{~}{\varphi }\varphi }^1Z_\varphi \mathrm{\Gamma }_{\stackrel{~}{\varphi }\varphi }|_{NP}=1.$$
(13)
Comparing these corrections with those appearing in the expression to $`\mathrm{\Gamma }_{\stackrel{~}{\varphi }\varphi \varphi }^{}`$, it is possible to see that the insertion of the composite operator $`\stackrel{~}{\varphi }\varphi `$ is equivalent to the “insertion” of the dashed-continuous-line leg in Fig. 2. By inspection of all possible Feynman diagrams, we conclude that this is valid to all orders in perturbation theory. Therefore the Feynman integrals involved in both $`\lambda ^{(1)}`$ and $`r`$ renormalizations are identical, and we can write:

$$\lambda ^{(1)}\mathrm{\Gamma }_{\stackrel{~}{\varphi }\varphi }=\mathrm{\Gamma }_{\stackrel{~}{\varphi }\varphi \varphi }^{}.$$
(14)
Consequently from Eq. (12), (13) and (14), it can be proved that,

$$\lambda _R^{(1)}=\lambda ^{(1)}Z_\varphi ^{1/2}Z_{\stackrel{~}{\varphi }\varphi }.$$
(15)
To describe how $`\lambda _R^{(1)}`$ flows under renormalization, it is necessary to define a Callan-Symanzik beta function. The dimensionless coupling constant corresponding to $`\lambda _R^{(1)}`$ is

$`g_R^{}`$ $`=`$ $`\left({\displaystyle \frac{\lambda _R^{(1)}}{D_R^2}}\right)\mu ^{2ϵ/2}`$ (16)
$`=`$ $`\left({\displaystyle \frac{\lambda ^{(1)}}{D_0^2}}\right)\mu ^{2ϵ/2}Z_\varphi ^{1/2}Z_{\stackrel{~}{\varphi }\varphi }Z_D^2,`$
where $`ϵ=4d`$, and $`D_R=D_0Z_D`$ is the renormalized diffusion coefficient. Then the beta function, which gives the differential renormalization group flow equation of $`g_R`$ is

$$\beta (g_R^{})=\mu \frac{g_R^{}}{\mu }=g_R^{}[2\frac{ϵ}{2}+\frac{1}{2}\gamma _\varphi +\gamma _{\stackrel{~}{\varphi }\varphi }2\gamma _D],$$
(17)
where $`\gamma _\varphi `$, $`\gamma _{\stackrel{~}{\varphi }\varphi }`$ and $`\gamma _D`$ are defined as

$$\gamma _\varphi =\mu \frac{lnZ_\varphi }{\mu },$$
(18)
$$\gamma _{\stackrel{~}{\varphi }\varphi }=\mu \frac{lnZ_{\stackrel{~}{\varphi }\varphi }}{\mu },$$
(19)
$$\gamma _D=\mu \frac{lnZ_D}{\mu }.$$
(20)
At the DP fixed point these gamma functions are related to the critical exponents by $`\gamma _\varphi ^{}=\frac{2\beta }{\nu _{}}d`$, $`\gamma _{\stackrel{~}{\varphi }\varphi }^{}=z\frac{1}{\nu _{}}`$ and $`\gamma _D^{}=z2`$, and thus

$$\beta (g_R^{})=\frac{1}{\nu _{}}(\beta \nu _{}1)g_R^{}+O(g_R^2).$$
(21)
Then, as we are approaching the infrared limit $`(k0)`$, $`g_R^{}`$ increases, being relevant to all orders in perturbation theory.

The renormalization group eigenvalue $`y_{\lambda ^{(1)}}`$, corresponding to $`\lambda ^{(1)}`$, is equal to

$$y_{\lambda ^{(1)}}=\frac{1}{\nu _{}}(1+\nu _{}\beta ).$$
(22)
Since this expression is always positive, the first term of the expansion in Eq. (9) is relevant to all orders of perturbation theory around $`d=d_c=4`$.

However, for $`d<d_c=4`$ the coupling constants are increasingly relevant in the expansion (9), and it is no longer meaningful to expand the exponential function. The system is driven away from the vicinity of the DP fixed point towards a non-trivial fixed point of order $`4d`$. To study this case, we apply a scaling analysis. Assume that position and time scale as $`[x]=k^1`$ and $`[t]=k^z`$. In the DP theory away from criticality, the non-zero mass term scales as $`[r]=k^{1/\nu _{}}`$. It couples to the operator $`\stackrel{~}{\varphi }\varphi `$, which has a scaling dimension $`x_{\stackrel{~}{\varphi }\varphi }`$, such that, $`[\stackrel{~}{\varphi }\varphi ]=k^{x_{\stackrel{~}{\varphi }\varphi }}`$. Using the fact that the action has to be dimensionless, we find that $`x_{\stackrel{~}{\varphi }\varphi }=d+z\frac{1}{\nu _{}}`$. Let us focus now on the scaling analysis of the non-Markovian term added to the DP action. First we consider the term corresponding to $`n=1`$. $`\lambda ^{(1)}`$ scales as $`[\lambda ^{(1)}]=k^{y_{\lambda ^{(1)}}}`$. Notice that the operators $`\stackrel{~}{\varphi }(x,t)\varphi (x,t)`$ and $`\varphi (x,t^{})`$ are evaluated at different space-time points, and consequently they are not correlated. So, each of them scales with its own scaling dimension, $`x_{\stackrel{~}{\varphi }\varphi }`$ and $`x_\varphi =\beta /\nu _{}`$, such that, $`[\stackrel{~}{\varphi }(x,t)\varphi (x,t)]=k^{x_{\stackrel{~}{\varphi }\varphi }}`$, $`[\varphi (x,t^{})]=k^{x_\varphi }`$. From the dimensionless nature of the action, we thus obtain $`y_\lambda =\frac{1}{\nu _{}}(\beta \nu _{}1)`$. Generalizing this for higher-order terms in the exponential expansion we find that the scaling exponent for each $`\lambda ^{(n)}`$ are,

$$y_{\lambda ^{(n)}}=\frac{1}{\nu _{}}\left(1+n(\nu _{}\beta )\right).$$
(23)
Since $`\nu _{}\beta `$ is always positive in $`d<4`$ dimensions, the terms of the expansion (9) are still increasingly relevant.

Furthermore, the scaling invariance of the non-Markovian term can only be established if the exponential function and its arguments are dimensionless. Therefore the two couplings $`w`$ and $`\lambda `$ scale separately with different scaling exponents, i.e., $`[\lambda ]=k^{y_\lambda }`$ and $`[w]=k^{y_w}`$. From Eq. (23) we obtain,

$$y_\lambda =\frac{1}{\nu _{}},y_w=\frac{\nu _{}\beta }{\nu _{}}.$$
(24)
Consequently the coupling constant $`w`$, is also relevant under renormalization.

## III Low first infection probability limit

In the previous section using scaling arguments, we showed that the non-Markovian term is relevant for $`d<d_c=4`$. We can now argue that a critical spreading process on the horizontal line in Fig. 1, is therefore driven away from the vicinities of the DP critical point to the points $`A`$ or $`B`$, where the probability of first infections is either very low or very high. In order to understand the dynamic behavior in these limits we study the asymptotic spreading properties keeping the second infection probabilities at the critical value of DP. We consider an initial state with a single infected site at the origin in a non-infected environment. For simplicity we will focus on the $`1+1`$-dimensional lattice. In this case there are not empty sites inside the epidemic clusters, which are composed of infected and immune sites only.

### III.1 Quasi-static approximation

Let us consider the limit of a very small first infection probability. In this regime, although first infections hardly ever happen, there still exists a surviving critical epidemic process as is shown in Fig. 4. The active domain is composed by re-infected sites among immune ones, and can be considered as bounded by almost rigid “walls” of empty sites. Let us define $`L(t)`$ to be the distance at time t between these two walls, such that there is at least one site active. If we keep the second infection rate critical, the dynamic epidemic process can be thought as an effective critical DP process evolving in a finite system of size $`L(t)`$.

A quasi-static approximation is based on the assumption that the time scale on which $`L(t)`$ grows is much larger than the correlation time of the DP process. Therefore, $`L(t)`$ can be assumed constant. It is known that in a finite-size system with a constant width $`L`$, the survival probability of a DP process decays exponentially as $`P(t)\mathrm{exp}(t/L^z)`$, where $`z=\frac{\nu _{}}{\nu _{}}`$. Then, in a quasi-static finite system the survival probability can be approximated as $`P_s(t)\mathrm{exp}(t/\xi _{})`$, i.e.

$$\frac{dP_s(t)}{dt}\frac{1}{\xi _{}}P_s(t)aP_s(t)L(t)^z,$$
(25)
where $`a`$ is a non-universal amplitude factor. We need now to know how $`L(t)`$ grows. Since the two boundaries are absorbing, we can apply the theory of DP with absorbing walls FHL01 . According to these results, the density of active sites next to the walls generated by the surviving clusters scales as

$$\rho _s(t)=bL(t)^{\beta _s/\nu _{}},$$
(26)
where $`b`$ is a non-universal amplitude factor which depends on the value of $`p_0`$. $`\beta _s0.73371(2)`$ is the surface critical exponent of DP in 1+1 dimensions. Clearly, $`dL/dt`$ is proportional to the frequency by which the system attempts to infect sites at the boundary and thus also to the density of active sites next to the boundary:

$$\frac{dL(t)}{dt}=p_0\rho _s(t)=p_0bL(t)^{\beta _s/\nu _{}}.$$
(27)
Therefore, the size of the domain grows as

$$L(t)=(p_0b(1+\beta _s/\nu _{})t)^{1/(1+\beta _s/\nu _{})}.$$
(28)
Inserting this result into Eq. (25) and solving the differential equation we obtain

$$\mathrm{ln}P_s(t)_0^tL(z)^z𝑑tp_0^\alpha t^{1\alpha },$$
(29)
where $`\alpha =\frac{\nu _{}}{\nu _{}+\beta _s}0.947167`$. Hence the survival decays asymptotically as a stretched exponential of the form

$$P_s(t)=P_s(0)\mathrm{exp}\left(Ap_0^\alpha t^{1\alpha }\right).$$
(30)
with $`A=\frac{a}{1\alpha }(\frac{\alpha }{bz})^\alpha `$. The above result implies that the average survival time $`T`$ is finite and scales as

$$Tp_0^{\alpha /(1\alpha )}.$$
(31)
It is interesting to compare these results with recent findings for random walkers between movable reflectors, where the survival probability was shown to decay as a power law with continuously varying exponents Dickman . We note that this result is not in contradiction with the present work since it would correspond to the limit $`\alpha 1`$.

### III.2 Numerical simulations

The lattice model considered is a 1+1 dimensional triangular lattice. We simulate a directed bond percolation process with a first infection rate $`p_0`$ different from a rate for subsequent infections or second infection rate $`p`$.

Let us first consider the limit of a very small first infection probability, where it is difficult for the process to infect sites at the boundary that have never been infected before (see Fig. 4). Bonds inside the infected region are open with the critical probability $`p=p_c=0.6447001(1)`$ Jensen96 , while bonds leading to healthy sites at the boundary are open with a different probability of first infection $`p_0`$. However, as $`p_0`$ is very small, conventional seed simulations are not suitable since most of the runs terminate after a very short time. For example, for $`p=0.1`$ most of the runs survive for less than $`100`$ time steps. Consequently we apply enrichment methods in order to circumvent this limitation. First of all, we apply a very simple enrichment method, in which we consider one lattice system. We keep the activity artificially alive at any time $`t`$. This is achieved by ignoring the updates which lead to the inactive state. Averaging over many realizations we measured the domain size $`L(t)`$ and the surface density of active sites. We find that the exponents predicted by the quasi-static approximation are in good agreement with the values obtained with the numerical simulations as shown in Fig. 5 and Fig. 6. The nonuniversal factor $`b`$ in Eq. (28) depends on the value of $`p_0`$. For, $`p_0=0.01`$ it takes the value $`b=1.65`$.

To determine the survival probability $`P_s(t)`$, we apply another enrichment method GrassWinners , which leads to a considerable improvement. The simulation starts with an ensemble of $`N=\mathrm{65\hspace{0.17em}536}`$ independent systems. Whenever the number of active runs becomes smaller than $`N/2`$ the ensemble is duplicated by creating identical copies of all the remaining active states. The new systems are labelled according to their ancestors. The survival probability $`P_s`$ then is reduced by a factor of $`n_f/n_i`$, where $`n_i`$ is the number of initial active systems and $`n_f`$ is the number of remaining active systems before the duplication. This process may be repeated as long as the ensemble has a sufficiently large number $`m`$ of independent ancestors at $`t=0`$. Using this method we were able to extend the temporal range of the simulation by four orders of magnitude up to $`t=10^6`$. Our results are shown in Fig. 7 for various values of $`p_0`$. As can be seen, the survival probability plotted in a double logarithmic representation is not a straight line, proving that it does not follow a power law. For the case $`p_0=0.01`$, the parameter $`a`$ in Eq. (25) takes the value $`a=1.00(2)`$.

### III.3 Lattice effects

The scaling arguments in Section III.1 are developed assuming the existence of a well defined mean value for the domain size $`L(t)`$. To justify this assumption, we define $`P(L,t)`$ to be the probability of having a domain of size $`L`$ at time $`t`$. Then the survival probability is $`P_s(t)=P(L,t)𝑑L`$. In this subsection we show that $`L(t)`$ is peaked around the mean value $`L(t)`$.

We begin by writing down a heuristic discrete master equation for the temporal evolution of $`P(L,t)`$:

$`P(L,t+1)`$ $`=`$ $`e^{aL^z}P(L,t)+bp_0(L1)^{\beta _s/\nu _{}}\times `$ (32)
$`P(L1,t)bp_0L^{\beta _s/\nu _{}}P(L,t).`$
The first term in the right hand side of Eq. (32) describes the change in $`P(L,t)`$ due to terminating runs according to Eq. (25) integrated over the time interval $`[t,t+1]`$. The second term corresponds to the probability of creating a domain of size $`L`$ from a domain of length ($`L1`$). The third term describes a loss term corresponding to $`LL+1`$. We neglect the second order contributions $`+b^2p_0^2(L1)^{2\beta _s/\nu _{}}P(L2,t)`$ and $`b^2p_0^2(L1)^{2\beta _s/\nu _{}}`$ which only affect the behaviour at small time and small $`L`$ regime.

It is straightforward to show that this master equation is consistent with Eq. (25), by summing over all $`L`$ and replacing $`L`$ by its average value $`L(t)`$. On doing this we obtain $`dP_s(t)/dt=aL(t)^zP_s(t)`$, that is the same equation as Eq. (25).

Consider Eq. (32) with $`p_0`$, $`a`$ and $`b`$ fixed at certain values. Then the differential equation which corresponds to the master equation, Eq. (32), can be written as

$`{\displaystyle \frac{P(L,t)}{t}}`$ $`=`$ $`{\displaystyle \frac{a}{L^z}}P(L,t)`$ (33)
$`p_0b{\displaystyle \frac{}{L}}\left[L^{\beta _s/\nu _{}}P(L,t)\right],L1.`$
The solution for Eq. (33) is

$`P(L,t)`$ $`=`$ $`e^{\frac{a}{p_0b(1z+\beta _s/\nu _{})}L^{1z+\beta _s/\nu _{}}}L^{\beta _s/\nu _{}}`$ (34)
$`\psi \left[L^{\beta _s/\nu _{}+1}p_0b(\beta _s/\nu _{}+1)t\right],`$
where $`L1`$. Consider the integrated distribution probability $`G(L,t)=_L^{\mathrm{}}P(L,t)𝑑L`$. Then $`G(L,t)`$ has the functional form

$`G(L,t)`$ $`=`$ $`e^{c\left[p_0b(\beta _s/\nu _{}+1)t\right]^\kappa }\times `$
$`\varphi \left[L^{\frac{\beta _s}{\nu _{}}+1}p_0b({\displaystyle \frac{\beta _s}{\nu _{}}}+1)t\right]\left(1+O({\displaystyle \frac{1}{t^{1\kappa }}})\right),`$
with $`c=a/(p_0b(1z+\beta _s/\nu _{}))`$ and $`\kappa =(1z+\beta _s/\nu _{})/(1+\beta _s/\nu _{})`$.

To observe the qualitative behaviour of $`P(L,t)`$ and $`G(L,t)`$, we carry out numerical integrations of Eq. (32). As initial conditions we choose the state in which there is only one site infected at time $`t=0`$. In Fig. 8, $`G(L,t)`$ is shown for several values of time $`t`$. It behaves like a propagating front with an overall exponential factor of time $`t`$. The fact that all the $`G(L,t)`$ curves intersect at the same point (see inset in Fig. 8), implies that there is a well defined peak and a mean value $`L(t)`$ for $`P(L,t)`$. The qualitative behaviour of $`G(L,t)`$ is independent of the particular values of $`a`$ and $`p_0`$. Consequently, to perform these calculations, we use values of $`a`$ and $`p_0`$ such that the effect of the exponential decay in Eq. (III.3) is not so pronounced.

We now compare the predictions of the master equation proposed in this subsection, with the results obtained from a simulation of the model. In Eq. (32), we fixed the value of $`p_0=0.01`$ and vary $`a`$. The survival probability calculated from direct numerical integration of Eq. (32), describes well the corresponding simulation result for $`p=0.01`$, in the case $`a=1`$ (see Fig. 9).

## IV High first infection probability limit

Let us finally consider the limit of a very high first infection probability $`p_01`$. Again we restrict to the 1+1-dimensional case, where the region of immune sites does not contain healthy sites inside. Because of the enhanced spreading probability at the boundaries, the domain grows rapidly by sudden avalanches of successive first infections (see Fig. 10). Since the avalanches at the left and the right boundary are expected to be uncorrelated, it suffices to study the propagation of one of the boundaries. As in the previous section, we propose simple scaling arguments in order to describe the growth of the infected/immune domain and the survival probability.

### IV.1 Independent avalanche approximation

In 1+1-dimensional directed bond percolation with $`p=p_c`$ and $`p_01`$ an avalanche is caused by a sequence of open bonds at the boundary. Therefore, the avalanche size $`\xi `$ is distributed exponentially as

$$P(\xi )=(1p_0)p_0^\xi e^{\xi /\overline{\xi }},$$
(36)
The quantity

$$\overline{\xi }=\frac{1}{\mathrm{ln}p_0}\frac{1}{1p_0}p_01.$$
(37)
is the average distance by which the avalanche advances the boundary in space. After each avalanche the process continues to evolve as an ordinary critical DP process inside the immune domain until it terminates or returns to the boundary where it releases a new avalanche. Thus the spreading behavior is mainly determined by the distribution of waiting times $`\tau `$ between the avalanches. We argue that the distribution of waiting times between avalanches is related to the problem of local persistence in DP HinrichsenKoduvely ; DuesseldorfPersistence  (for a recent review on persistence see e.g. Satya99 ). The local persistence probability $`R(t)`$ is defined as the probability that a randomly selected site in an ordinary critical DP process starting from a homogeneous initial state has not been reactivated until time $`t`$. It was shown that this quantity decays algebraically as $`R(t)t^\mathrm{\Theta }`$, where $`\mathrm{\Theta }=1.50(1)`$ is the so-called local persistence exponent HinrichsenKoduvely . In the present problem the situation is similar: Each avalanche creates locally a quasi-homogeneous state. The process then evolves as an ordinary critical DP process inside the infected/immune region until the boundary is revisited for the first time in order to release a new avalanche. However, unlike persistence studies in 1+1 dimensions, where a persistent site can be activated independently from the left and from the right, the boundary sites in the present problem can be infected only from one side. Hence the probability that the next avalanche has not yet been released decays as $`\tau ^{\mathrm{\Theta }/2}`$. Thus we conjecture that the waiting times between avalanches are distributed algebraically as

$$P(\tau )\tau ^{1\mathrm{\Theta }/2}.$$
(38)
Next, we argue that these waiting times may be interpreted as directed Lévy flights in time LevyGeneral . After each flight the domain size $`L(t)`$ grows on average by the mean avalanche size $`\overline{\xi }`$. This type of growth may be described by the equation

$$D_t^{\mathrm{\Theta }/2}L(t)\overline{\xi },$$
(39)
where $`D_t^{\mathrm{\Theta }/2}`$ is a fractional derivative defined through its action in Fourier space $`D_t^{\mathrm{\Theta }/2}e^{i\omega t}=(i\omega )^{\mathrm{\Theta }/2}e^{i\omega t}`$. Simple dimensional analysis leads to the result

$$L(t)\overline{\xi }t^{\mathrm{\Theta }/2}.$$
(40)
In order to compute the survival probability, we note that during the waiting time the outermost active site departs from the boundary with an average distance $`\mathrm{}(\tau )\tau ^{1/z}`$. Obviously, the process can only terminate if this distance is of the same order as the size of the infected/immune domain $`L(t)`$. Therefore, we expect the distribution of waiting times between avalanches to be cut off by a maximal waiting time

$$\tau _{\text{max}}L^z(t)\overline{\xi }^zt^{z\mathrm{\Theta }/2}.$$
(41)
Consequently, the probability that the process terminates between two avalanches is given by

$$P_0=\frac{_{\tau _{\text{max}}}^{\mathrm{}}𝑑\tau P(\tau )}{_1^{\mathrm{}}𝑑\tau P(\tau )}\tau _{\text{max}}^{\mathrm{\Theta }/2}$$
(42)
On the other hand, the cutoff due to terminating runs implies that the average waiting time is finite and scales as

$$\overline{\tau }=\frac{_1^{\tau _{\text{max}}}𝑑\tau \tau P(\tau )}{_1^{\tau _{\text{max}}}𝑑\tau P(\tau )}=\tau _{\text{max}}^{1\mathrm{\Theta }/2}$$
(43)
Therefore, the average loss of the survival probability per unit time is given by $`P_0/\overline{\tau }1/\tau _{\text{max}}`$, i.e.

$$\frac{dP_s(t)}{dt}P_s(t)/\tau _{\text{max}}P_s(t)\overline{\xi }^zt^{z\mathrm{\Theta }/2}.$$
(44)
Solving this equation, the asymptotic behavior of the survival probability is not a stretched exponential and is given by a

$$P_s(t)=\frac{1}{\stackrel{~}{N}}\mathrm{exp}\left(+\stackrel{~}{A}\overline{\xi }^zt^{1z\mathrm{\Theta }/2}\right),$$
(45)
where $`\stackrel{~}{N}`$ and $`\stackrel{~}{A}`$ are unknown constants. In contrast to the previous case in Eq. (30) the exponent $`1z\mathrm{\Theta }/20.185`$ is negative. Consequently the survival probability tends to a constant $`P_s(\mathrm{})>0`$, meaning that a finite fraction of runs survives for infinitely long time. For very large $`p_0`$ this constant is expected to scale as

$$1P_s(\mathrm{})(1p_0)^{z\mathrm{\Theta }/2}.$$
(46)
### IV.2 Numerical results

In order to verify these results numerically, we developed an especially optimized Monte Carlo algorithm for large values of $`p_0`$. In order to compute the survival time for a run and to see how the boundaries advance for a given realization of open and closed bonds, it is in most cases not necessary to construct the entire cluster. Rather it suffices to construct its branches next to the boundaries. For example, in Fig. 10 large parts of the cluster are irrelevant for the advancement of the boundaries and the survival time. For this reason our algorithm first constructs the branches of the cluster next to the boundaries. If these branches terminate before the simulation time is reached, we recursively constructs the omitted branches in the interior of the clusters. If however they do not terminate, then the recursive construction is not performed. Using this technique we could extend the simulation time by two decades to $`10^4`$ time steps, two decades less than in the previous case of low $`p_0`$.

Our numerical results are shown in Fig. 11. The left panel shows the survival probability as a function of time in a directed bond percolation process with $`p_0=0.85,0.90,0.95,`$ and $`0.98`$. The positive curvature of the lines indicates that there is no power-law scaling. Although the upward curvature is in agreement with the expected result (45), the simulation time is not large enough to confirm the behaviour of Eq. (45) quantitatively, mainly because the constant $`\stackrel{~}{A}`$ is not known. In order to substantiate the assumptions made in the previous subsection, we first measured the growing domain size $`L(t)`$. As shown in Fig. 11b, the measured slope seems to tend to the predicted slope $`\mathrm{\Theta }/2=0.75(1)`$. Moreover, we estimated the logarithmic derivative of the survival probability, in which the unknown prefactor $`\stackrel{~}{A}`$ drops out (see right panel of Fig. 11). Although this data set is quite noisy, we observe a rather clean power law. The estimated slope $`1.12(10)`$ is consistent with the theoretical prediction $`t^{z\mathrm{\Theta }/2}=t^{1.185}`$, supporting the assumptions which led to the result (45).

## V Conclusions

We analyzed the effects of immunization as a small perturbation on the DP model and studied in detail the scaling behaviour of the theory around the DP critical point. We derived by an alternative method the field theoretic action for the model. A non-Markovian term is added to the DP action because of the presence of immunization. In the field theory, the probability for subsequent infections, which is different from the first infection probability, is a function of position and time. Nevertheless, we assumed in the lattice model that the susceptibility for spreading changes only after the first infection, remaining constant thereafter. This assumption does not changes the final results.

The phase diagram (cf. Fig. 1) comprises two phase transition lines, namely, a horizontal line, where the process in reinfected regions shows the critical behavior of DP, and (in more than one spatial dimension) a curved transition line, where the critical behavior corresponds to the general epidemic process studied in GCR97 . Both lines meet at the DP critical point $`p=p_0=p_c`$.

The non-Markovian term turns out to be relevant for $`d<6`$. We considered the RG flows and argued that any point in the neighbourhood of the DP critical point will be driven away from it. In particular, we focused our study on the horizontal phase transition line with a critical reinfection probability $`p=p_c`$. The system along this line, is driven away from the DP critical point as soon as the immunization effects are turned on as a small perturbation. The asymptotic behavior is determined by the limits of very low and high first infection probability (close to points $`A`$ and $`B`$ in Fig. 1).

We proposed simple scaling arguments for the behaviour of the survival probability in both limits. We related the corresponding exponents with those exponents of the critical ordinary DP theory and DP near a wall. The survival probability is found to obey a stretched exponential behaviour in the low first infection probability limit, and it decays to a constant in the high first infection probability limit, taking the form

$$P_s(t)\{\begin{array}{cc}\mathrm{exp}(at^{+0.0528})\hfill & \text{ for }p_00\hfill \\ \mathrm{exp}(+bt^{0.185})\hfill & \text{ for }p_01\hfill \end{array}$$
(47)
The numerical simulations in $`1+1`$ dimensions support these theoretical predictions, ruling out the possibility of asymptotic power-law scaling in both limits. The question as to whether the stretched exponential decay of the survival probability persists in higher dimensions is still open.

Finally, we comment on a possible connection of these results with the problem of infinitely many absorbing states. A typical example of such systems is the pair contact process $`2A3A`$, $`2A\mathrm{}`$, in which solitary particles are not allowed to diffuse Jensen93 . In Ref. Munoz96 ; LopezMunoz97 ; Munoz97 ; MGD98  this model was studied with an effective Langevin equation and it was inferred that the survival probability decays as a power law at the transition point, with continuously varying exponents. However, the field theoretical action studied in Munoz96 ; Munoz97 ; MGD98  is identical to the one studied in this paper. Hence, we conclude that, according to the analysis presented in this paper, seed simulations of the pair contact process, at least in $`1+1`$ dimensions, should show an stretched exponential behaviour and not a power low decay at the transition point.

###### Acknowledgements.

We want to thank Peter Grassberger for communication of unpublished numerical calculations. A.J.D. is greatly indebted to John Cardy for enlightening discussions and suggestions, and for his most careful reading of the manuscript. A.J.D. would also like to thank to R. Rajesh for useful comments. Part of the present computations have been carried out on the 128-node Alpha Linux Cluster Engine ALiCE at Wuppertal University. We thank N. Eicker, Th. Lippert and B. Orth for their assistance. A.J.D. was supported by CONICET (Argentina), the British Council-Fundación Antorchas (Argentina), and the ORS Award Scheme (UK).


text:
# Quantum walks on cycles

## I  Introduction

Celebrated results of Shor Shor (1994) and Grover Grover (1977) started a quest for algorithms based on quantum mechanics which can surpass corresponding classical algorithms. As many classical algorithms employ properties of random walks on graphs, several groups of researchers have begun to study properties of quantum analogues of classical random walks (for a survey see Kempe (a)). Two models of quantum random walks have been proposed by Aharonov et al. Aharonov et al.  and Farhi and Gutmann Farhi ; in the paper we consider only a discrete quantum walk as defined in Aharonov et al. . The behavior of quantum random walks have been shown to differ greatly from their classical counterparts. Ambainis et al. Ambainis et al.  proved that the spreading time in the quantum walk on line scales linearly with the number of steps, while Aharonov et al. Aharonov et al.  showed that the mixing time for a walk on a cycle grows linearly with the cycle length. Even more spectacular exponential speed up was discovered by Kempe et al. Kempe (b) who studied a quantum random walk on hypercube; this result led to the construction of the first quantum algorithm based on a random walk by Shenvi et al. Shenvi et al. . The fact that in quantum walks on graphs the probability function spreads out much faster than in classical case is not the only factor which can be explored in designing new quantum algorithms. In Aharonov et al.  the authors remarked that ‘one may try to use quantum walks which converge to limiting distribution which are different than those of the corresponding classical walks’. In this note we follow this suggestion and study the limiting distribution of a random walk on a cycle; it turns out that it depends on the length of the cycle in a somewhat surprising way. We also give an example of a quantum walk in which the distance from the initial state to the uniform distribution is larger than the distance between the uniform distribution and the initial state. Finally, we remark that recently Travaglione and Milburn Travaglione and Milburn (2002) and Dür et al. Dür et al. (2002) have proposed a scheme of experiment which realizes a quantum walk on a cycle.

## II  Model

We study a quantum random walk on a cycle with $`d`$ nodes. In the model of such a walk proposed in Aharonov et al.  nodes are represented by vectors $`|v`$, $`v=0,1,\mathrm{},d1`$, which form an orthonormal basis of the Hilbert space $`H_V`$. An auxiliary two-dimensional Hilbert space $`H_A`$ (coin space) is spanned by vectors $`|s`$, $`s=0,1`$. The initial state of the walk is a normalized vector

$$|\mathrm{\Psi }_0=\underset{s,v}{}\alpha _{sv}|s,v=\underset{s,v}{}\alpha _{sv}|s|v$$
(1)
from the tensor product $`H=H_AH_V`$. In a single step of the walk the state changes according to the equation

$$|\mathrm{\Psi }_{n+1}=U|\mathrm{\Psi }_n,$$
(2)
where the operation $`U=S(HI)`$ first applies the Hadamard gate operator $`H=_{s,s^{}}(1)^{ss^{}}|ss^{}|`$ to the vector from $`H_A`$, and then shifts the state by the operator

$$S=\underset{s,v}{}|ss||v+2s1(\mathrm{mod}d)v|.$$
(3)
The probability distribution on the nodes of the cycle after the first $`n`$ steps of the walk is given by

$$p_n(v)=\underset{s}{}|s,v|\mathrm{\Psi }_n|^2.$$
(4)
However, as was observed by Aharonov at el. Aharonov et al. , for a fixed $`v`$, the probability $`p_n(v)`$ is ‘quasi-periodic’ as a function of $`n`$ and thus, typically, it does not converge to a limit. Thus, instead of $`p_n(v)`$ the authors of Aharonov et al.  considered

$$\overline{p}_n(v)=\frac{1}{n}\underset{i=1}{\overset{n}{}}p_n(v),$$
(5)
and proved that for any initial state $`|\mathrm{\Psi }_0`$ and every node $`v=0,1,\mathrm{},d1`$, the sequence $`\overline{p}_n(v)`$ converges to the limiting distribution

$$\pi (v)=\underset{a,a^{}}{}\underset{s}{}\mathrm{\Gamma }_{aa^{}}|\varphi _a|\mathrm{\Psi }_0\mathrm{\Psi }_0|\varphi _a^{}s,v|\varphi _a\varphi _a^{}|s,v,$$
(6)
where $`c_a`$ are the eigenvalues of $`U`$, $`|\varphi _a`$ stand for the eigenvectors of $`U`$, and

$$\mathrm{\Gamma }_{aa^{}}=\{\begin{array}{cc}1\hfill & \mathrm{if}c_a=c_a^{}\hfill \\ 0\hfill & \mathrm{if}c_ac_a^{}.\hfill \end{array}$$
(7)
In the case of the Hadamard walk on cycle, for $`j=0,1`$ and $`k=0,1,\mathrm{},d1`$, we get

$$c_{jk}=\frac{1}{\sqrt{2}}\left((1)^k\sqrt{1+\mathrm{cos}^2(2\pi j/d)}i\mathrm{sin}(2\pi j/d)\right),$$
(8)
and

$$|\varphi _{jk}=(a_{jk}|0+a_{jk}b_{jk}|1)\underset{v}{}\omega _d^{jv}|v,$$
(9)
where $`\omega _d=e^{2\pi i/d}`$,

$`a_{jk}`$ $`=1/\sqrt{d(1+|b_{jk}|^2)},`$ (10)
$`b_{jk}`$ $`=\omega _d^j\left((1)^k\sqrt{1+\mathrm{cos}^2(2\pi j/d)}\mathrm{cos}(2\pi j/d)\right).`$ (11)
Note that $`c_{j0}c_{j^{}1}`$ so $`\mathrm{\Gamma }_{jk,j^{}k^{}}=\gamma _{jj^{}}\delta _{kk^{}}`$, where $`\delta _{kk^{}}`$ is the Kronecker delta. Thus, (6) becomes

$$\pi (v)=\underset{j,j^{}}{}\underset{k}{}\gamma _{jj}|\varphi _{jk}|\mathrm{\Psi }_0\mathrm{\Psi }_0|\varphi _{j^{}k}A_{jj^{}k}\omega _d^{v(jj^{})},$$
(12)
where $`A_{jj^{}k}=a_{jk}a_{j^{}k}(1+b_{jk}b_{j^{}k}^{})`$. If $`d`$ is odd, then all eigenvalues are distinct, $`\gamma _{jj^{}}=\delta _{jj^{}}`$, and $`A_{jjk}=1/d`$; consequently, $`\pi (v)=(1/d)_j_k|\varphi _{jk}|\mathrm{\Psi }_0|^2=1/d`$. It comes as no surprise, since as was proved in Aharonov et al. Aharonov et al.  the limiting distribution $`\pi `$ is always uniform in a non-degenerate case. Thus we concentrate on more interesting case of even $`d`$. Then, because of symmetries $`c_{d/2j,k}=c_{j,k}`$ and $`c_{d/2+j,k}=c_{jk}^{}=c_{dj,k}`$, the coefficient $`\gamma _{jj^{}}`$ does not vanish when one of the following conditions hold:

1. $`j=j^{}`$;
2. $`j=0`$, $`j^{}=d/2`$;
3. $`j=t`$, $`j^{}=d/2t`$, for $`t=1,2,\mathrm{},t_{\mathrm{max}}`$;
4. $`j=dt`$, $`j^{}=d/2+t`$, for $`t=1,2,\mathrm{},t_{\mathrm{max}}`$;

where here and below $`t_{\mathrm{max}}=(d2)/4`$.

## III Result and discussion

We study in detail the limiting distribution $`\pi `$ for a quantum walk for which starts with a state in which with probability one a particle is at a node $`v_0`$; more specifically we set

$$|\mathrm{\Psi }_0=\frac{1}{\sqrt{2}}\left(|0,v_0+i|1,v_0\right).$$
(13)
Then $`\varphi _{jk}|\mathrm{\Psi }_0=g_{jk}\omega _d^{v_0j}`$, where $`g_{jk}=a_{jk}(1+ib_{jk}^{})/\sqrt{2}`$, and

$$\pi (v)=\underset{j,j^{}}{}\underset{k}{}\gamma _{jj^{}}f(j,j^{},k)\omega _d^{(vv_0)(jj^{})},$$
(14)
with $`f(j,j^{},k)=g_{jk}g_{j^{}k}^{}A_{jj^{}k}`$. Since $`\gamma _{jj^{}}=0`$ except of the four cases (i)–(iv) described above, we get

$$\pi (v)=F+(1)^\mathrm{\Delta }2\mathrm{R}\mathrm{e}(F_0)+(1)^\mathrm{\Delta }\mathrm{Re}\left(\underset{i=1}{\overset{t_{\mathrm{max}}}{}}F_t\right),$$
(15)
where by

$$\mathrm{\Delta }=\mathrm{\Delta }(v)=\mathrm{min}\{|vv_0|,d|vv_0|\}$$
(16)
denotes the distance between nodes $`v_0`$ and $`v`$, and

$`F`$ $`={\displaystyle \underset{j,k}{}}f(j,j,k),`$ (17)
$`F_0`$ $`={\displaystyle \underset{k}{}}f(0,d/2,k),`$
$`F_t`$ $`={\displaystyle \underset{k}{}}f(t,d/2t,k)+f(d/2+t,dt,k).`$
After some elementary but not very exciting calculations (17) reduces to

$`F`$ $`={\displaystyle \frac{1}{d}},`$ (18)
$`F_0`$ $`=0,`$
$`F_t`$ $`={\displaystyle \frac{2}{d^2}}{\displaystyle \frac{\mathrm{sin}^2(2\pi t/d)}{1+\mathrm{cos}^2(2\pi t/d)}}.`$
We remark that $`F_t=0`$ whenever $`t_{\mathrm{max}}=0`$, i.e., for $`d=2`$ and $`d=4`$. Hence the limiting distribution for even cycles of sizes two and four are uniform, which has also been observed by Travaglione and Milburn Travaglione and Milburn (2002), who analyzed a quantum walk on cycle of length four step by step. However, for $`d6`$ we have

$$\pi (v)=(1+\mathrm{\Pi }(v))/d,$$
(19)
with the ‘correction’ term $`\mathrm{\Pi }(v)=(1)^\mathrm{\Delta }4S/d`$, where

$$S=\underset{t=1}{\overset{t_{\mathrm{max}}}{}}\mathrm{cos}(4\pi t\mathrm{\Delta }/d)\left(\frac{2}{1+\mathrm{cos}^2(2\pi t/d)}1\right).$$
(20)
Setting $`z=32\sqrt{2}0.17157\mathrm{}`$, one can write (20) as

$`S=\left({\displaystyle \frac{8z}{1z^2}}1\right){\displaystyle \underset{t=1}{\overset{t_{\mathrm{max}}}{}}}\mathrm{cos}(4\pi \mathrm{\Delta }t/d)`$ $`+{\displaystyle \frac{8z}{1z^2}}{\displaystyle \underset{m=1}{\overset{\mathrm{}}{}}}(z)^m{\displaystyle \underset{t=1}{\overset{t_{\mathrm{max}}}{}}}\mathrm{cos}(4\pi (\mathrm{\Delta }+m)t/d)`$ (21)
$`+{\displaystyle \frac{8z}{1z^2}}{\displaystyle \underset{m=1}{\overset{\mathrm{}}{}}}(z)^m{\displaystyle \underset{t=1}{\overset{t_{\mathrm{max}}}{}}}\mathrm{cos}(4\pi (\mathrm{\Delta }m)t/d),`$
which, in turn, transforms to

$`S={\displaystyle \frac{d}{4}}(\delta _{\mathrm{\Delta },0}+\delta _{\mathrm{\Delta },d/2})+`$ $`(1)^\mathrm{\Delta }{\displaystyle \frac{d}{4}}{\displaystyle \frac{8z}{1z^2}}{\displaystyle \frac{z^\mathrm{\Delta }(z)^{d/2}+z^\mathrm{\Delta }}{1(z)^{d/2}}}`$ (22)
$`+{\displaystyle \frac{1+(1)^\mathrm{\Delta }\xi }{2}}{\displaystyle \frac{4z}{(1z)^2}}{\displaystyle \frac{4z(1)^\mathrm{\Delta }\xi }{(1+z)^2}},`$
where $`\xi =\left(1+(1)^{d/2}\right)/2`$, i.e., $`\xi =1`$ when $`d/2`$ is even, and $`\xi =0`$ if $`d/2`$ is odd. Thus, we arrive at

$`\mathrm{\Pi }(v)=(1)^{\mathrm{\Delta }+1}(\delta _{\mathrm{\Delta },0}+\delta _{\mathrm{\Delta },d/2})+`$ $`{\displaystyle \frac{8z}{1z^2}}{\displaystyle \frac{z^\mathrm{\Delta }(z)^{d/2}+z^\mathrm{\Delta }}{1(z)^{d/2}}}`$ (23)
$`+{\displaystyle \frac{4}{d}}\left({\displaystyle \frac{\xi +(1)^\mathrm{\Delta }}{2}}{\displaystyle \frac{4(1)^\mathrm{\Delta }z}{(1z)^2}}{\displaystyle \frac{4z\xi }{(1+z)^2}}\right).`$
In Figure 1 we pictured the resulting limiting distributions for $`d=24`$ and $`d=26`$. It is easy to see they are almost uniform except for the nodes which lie next to the initially populated node $`v_0`$ and the opposite node $`\widehat{v}_0=v_0+d/2`$ $`(\mathrm{mod}d)`$. Indeed, as $`d\mathrm{}`$, the last term of (23) vanishes and

$$\overline{\mathrm{\Pi }}(v)=\underset{d\mathrm{}}{lim}\mathrm{\Pi }(v)=\eta (\mathrm{\Delta })(1)^\xi \eta (\mathrm{\Delta }^{}),$$
(24)
where

$$\eta (x)=\frac{8z^{1+x}}{1z^2}\delta _{x0},$$
(25)
and $`\mathrm{\Delta }^{}`$ is the distance between nodes $`v`$ and $`\widehat{v}_0`$. The equation (24) shows that the correction term $`\overline{\mathrm{\Pi }}(v)`$ is significant only for $`v`$ which lie close to either $`v_0`$ or $`\widehat{v}_0`$ and decreases exponentially with the distance between $`v`$ and $`v_0`$ and $`\widehat{v}_0`$. Note that the shapes of the cusps near $`v_0`$ and $`\widehat{V}_0`$ does not depend very much on $`d`$ for large $`d`$ (except of the scaling factor $`1/d`$). However, if $`d/2`$ is odd then the limiting distribution $`\pi (v)`$ has a minimum at $`\widehat{v}_0`$, while if $`d/2`$ is even the distribution has a peak $`\widehat{v}_0`$, virtually identical with that which appear at $`v_0`$ (see Figures 1 and 2). The fact that such a local behaviour of the limiting distribution $`\pi `$ depends so strongly on ‘global’ properties of space, as the parity of $`d/2`$, is somewhat surprising. We hope that this and/or analogous phenomena can be used in constructing efficient quantum algorithms.

It is easy to see that, despite of a narrow cusps near nodes $`v_0`$ and $`\widehat{v}_0`$, the total variation distance $`d_{TV}`$ between the limiting distribution given by (19) and (23) and the uniform distribution tends to 0 as $`d\mathrm{}`$. However it is not hard as well to construct ‘highly non-uniform’ distributions which remain invariant during the walk: one can simply take as the initial state a superposition of two degenerated eigenvectors (the distribution corresponding to a single eigenvector is always a uniform one). An example of such an invariant nonuniform probability distribution is presented on Figure 3.

With a slightly more work one can construct the initial states which indicate differences between quantum and classical walks on cycles even more distinctively. One such example is shown on Figure 4. Here the total variation distance from the uniform distribution, defined as

$$\frac{1}{2}\underset{v=0}{\overset{d1}{}}\left|p(v)\frac{1}{d}\right|$$
(26)
is equal $`0.046`$ for the initial state, while for the limiting distribution it grows to $`0.204`$.

## IV  Conclusion

We study the properties of a Hadamard quantum walk on a cycle with $`d`$ nodes. In the case of a walk starting from a single node we give an explicit formula for the limiting distribution and show that it is very sensitive to the arithmetic properties $`d`$. We hope that this or an analogous mode of behaviour can be used in construction of efficient quantum algorithms.

Furthermore, we present an example of a quantum walk on a cycle, for which the total variation distance between the initial distribution and the uniform distribution is much smaller than the distance between limiting distribution and the uniform one.

###### Acknowledgements.

We wish to thank the State Committee for Scientific Research (KBN) for its support: M.B. and T.Ł. were supported by grant 2 P03A 016 23; A.G. and A.W. by grant 0 T00A 003 23.