text:
# Statistics of the one-dimensional Riemann walk

## 1 Introduction

The Lévy flight is a random walk in continuous space whose step size distribution has a power law tail and is therefore sometimes called a ”Lévy  distribution”. The ubiquity of such distributions has been emphasized by many authors, and is a consequence of the power law tail being invariant under convolution. Many interesting instances of the occurrence of Lévy distributions are given by Tsallis  and Tsallis et al. . These range from applications in physics (superdiffusion, chaotic fluid flow) and engineering (leaking taps) through studies of the physiology of heart activity, all the way to descriptions  of fluctuations of financial markets.

A one-dimensional lattice version of the Lévy flight may be constructed as follows. Let a random walk consist of independent steps, and let the probability $`p(\mathrm{})`$ for a displacement of $`\mathrm{}`$ lattice units in a single step be given by $`p(0)=0`$ and

$$p(\mathrm{})=A|\mathrm{}|^{1\alpha }(\mathrm{}=\pm 1,\pm 2,\mathrm{})$$
(1.1)
Here $`\alpha >0`$ is the Lévy exponent. Normalization of $`p`$ implies that $`A^1=2\zeta (1+\alpha )`$ where $`\zeta `$ is the Riemann zeta function. This random walk was first studied by Gillis and Weiss  in 1970. It is called the Riemann walk by Hughes (Ref. , p. 154) and we will conform to that terminology. More generally we call of Riemann type any one-dimensional lattice walk whose $`p(\mathrm{})`$ is asymptotically proportional to $`|\mathrm{}|^{1\alpha }`$ when $`|\mathrm{}|\mathrm{}`$.

Riemann type walks were reviewed in detail by Hughes . Of particular interest is the exponent regime $`0<\alpha 2`$, where these walks have a mean square displacement per step, $`\mathrm{}^2`$, which is infinite. There then exists, at least for certain global walk features, a correspondence between simple random walk on a $`d`$-dimensional lattice and one-dimensional Riemann type walks of exponent $`\alpha =2/d`$. In some ways the fraction $`\frac{2}{\alpha }`$ acts as the walk’s effective dimensionality. But whereas analytical results for noninteger dimension $`d`$ cannot be checked by computer simulations, the full continuum of $`\alpha `$ values is accessible to Monte Carlo studies.

Much interest has centered around the following question. Let there be a $`t`$ step Riemann walk. Then what are the statistical properties of its support $`𝒮(t)`$, i.e., of the set of sites that the walk has visited? There appears immediately an important difference between the exponent regimes $`0<\alpha <1`$ and $`1<\alpha <2`$. In the former regime the Riemann walk is transient  and it is easy to show (see Sec. 2.5) that $`𝒮(\mathrm{})`$ is a set of fractal dimension $`d_𝒮=\alpha `$. In the latter case the Riemann walk is recurrent , $`𝒮(\mathrm{})`$ coincides with the full one-dimensional lattice, and $`d_𝒮=1`$. The existing literature deals with the different question of finding the properties of $`𝒮(t)`$ for asymptotically $`t`$; the results reflect, nevertheless, the same distinction between $`0<\alpha <1`$ and $`1<\alpha <2`$. The borderline case $`\alpha =1`$ is more subtle.

Gillis and Weiss  study the number $`S(t)`$ of distinct sites in the support. They find, among other results, that for $`t\mathrm{}`$ the average of this random variable behaves  as $`\overline{S(t)}t`$ for $`0<\alpha <1`$ and as $`\overline{S(t)}t^{1/\alpha }`$ for $`1<\alpha <2`$, where $``$ indicates asymptotic proportionality. For $`\alpha =1`$ and $`\alpha =2`$ power laws with logarithmic correction factors appear . For $`\alpha >2`$ the result $`\overline{S(t)}t^{1/2}`$ is identical to that for the simple random walk in $`d=1`$.

A recent extension of this work is due to Berkolaiko et al. . Pursuing a question initially asked for the case of the simple random walk by Larralde et al. , these authors investigate the number $`S_N(t)`$ of distinct sites visited by $`N`$ independent $`t`$ step Riemann type walks all starting on the same lattice site. Again power laws appear, both for $`t\mathrm{}`$ at fixed $`N`$ and for $`N\mathrm{}`$ at fixed $`t`$.

The present work extends the investigations of Gillis and Weiss into a different direction. We limit ourselves to the Riemann walk defined by Eq. (1.1), with $`\alpha `$ in the regime of greatest interest, that is, $`0<\alpha <2`$. Our results may be summarized under three headings.

1. Variance $`\overline{\mathrm{\Delta }S^2(t)}`$. For any quantity $`X(t)`$ we will denote its instantaneous deviation from average by $`\mathrm{\Delta }X(t)X(t)\overline{X(t)}`$. Traditionally in this field the calculation of the average number $`\overline{S(t)}`$ of distinct sites visited has been followed by a calculation of the variance $`\overline{\mathrm{\Delta }S^2(t)}`$ of that number. Thus, for the simple random walk $`\overline{S(t)}`$ was first calculated by Dvoretzky and Erdös  in 1951, and $`\overline{\mathrm{\Delta }S^2(t)}`$ by Jain and Pruitt  in 1970. For the one-dimensional Riemann walk the present work supplements the 1970 results due to Gillis and Weiss  for $`\overline{S(t)}`$ by the corresponding ones for the variances $`\overline{\mathrm{\Delta }S^2(t)}`$ in the regime $`0<\alpha <2`$.

2. Variables other than $`S(t)`$. We use the powerful generating function method (GFM), which was introduced into the field of random walks by Montroll  and Montroll and Weiss . Overviews of this method are given by Weiss  and by Hughes . The first calculation of a variance by the GFM, viz. that of $`S(t)`$ for the simple random walk, is due to Torney  in 1986.

In 1994 Coutinho et al.  performed Monte Carlo simulations of, among other things, the number of unvisited islands enclosed by the support of the $`t`$ step simple random walk in two dimensions. This led Caser and Hilhorst  to analytically determine the asymptotic behavior of the average number of islands. Subsequently Van Wijland et al.  developed a compact GFM based analytical scheme for calculating simultaneously the averages, variances, and correlations of a large class of observables characteristic of the support, generically denoted by the symbol $`M(t)`$. In $`d=2`$ this class includes also the total boundary length of the support, and in $`d=3`$ its surface area and Euler index.

Here we bring this scheme to bear on the one-dimensional Riemann walk. The support of this walk consists of alternating sequences of visited and unvisited sites. Among the most prominent members of the class of observables $`M(t)`$ is, next to $`S(t)`$, the number of visited sequences, that we will denote by $`I(t)`$. Islands in $`d=1`$ just are unvisited sequences enclosed by the support, of which there are $`I(t)1`$; the support furthermore has $`2I(t)`$ boundary sites (= visited sites adjacent to an unvisited one). Table I summarizes our results for the asymptotic laws of the averages, variances, and correlations involving $`S(t)`$ and $`I(t)`$. Beyond their intrinsic interest these laws may serve in heuristic arguments in reaction–diffusion processes, e.g., to estimate the trapping probability of an atom that diffuses in a random absorbing environment, or the effective reaction rate between two diffusing species. We defer further comments to Sec. 7.

3. Universality of fluctuations. The deviations from average $`\mathrm{\Delta }S(t)`$,$`\mathrm{\Delta }I(t)`$, $`\mathrm{}`$,$`\mathrm{\Delta }M(t)`$,$`\mathrm{}`$ are randomly time-dependent variables that one would a priori expect to exhibit some degree of correlation. One calls these fluctuations universal – by lack of a better name – when in the limit $`t\mathrm{}`$ all $`\mathrm{\Delta }M(t)`$ are asymptotically equal (up to a proportionality constant) to a single $`M`$ independent stochastic process. For the simple random walk universality was shown to hold in dimensions $`d=2`$  and $`d=3`$ , and not to hold in $`d=4,5,\mathrm{}`$. For the $`d=1`$ Riemann walk we find that universality holds in the exponent regime $`\frac{2}{3}\alpha <2`$, but not for $`0<\alpha <\frac{2}{3}`$. A novelty with respect to the case of the simple random walk is that for $`1<\alpha <2`$ not a single, but two $`M`$ independent processes are needed to describe the universal fluctuations: one applies to bulk and the other to surface observables. The precise statements are given in Sec. 6.

This article is set up as follows. Sec. 2 describes those elements of our analysis that are common to the full exponent interval $`0<\alpha <2`$. Secs. 3 and 4 deal more in particular with the exponent regimes $`0<\alpha <1`$ and $`1<\alpha <2`$, respectively, and derive the asymptotic behavior of averages, variances, and correlations. In Sec. 5 we do the same for the exceptional values $`\alpha =\frac{2}{3}`$ and $`\alpha =1`$. In Sec. 6 we discuss the universality properties. In Sec. 7 we provide some additional interpretation of our results and conclude.

## 2 Observables, averages, and correlations

### 2.1 Observables $`M(t)`$

Our analysis is based on first writing quantities of interest in terms of the field $`m(x,t)`$ of ”complementary occupation numbers” defined by $`m(x,t)=1`$ if site $`x`$ has not yet been visited at time $`t`$, and $`m(x,t)=0`$ otherwise. The expressions of $`S`$ and $`I`$ in terms of $`m`$ are

$`S(t)`$ $`=`$ $`{\displaystyle \underset{x=\mathrm{}}{\overset{\mathrm{}}{}}}[1m(x,t)]`$ (2.1)
$`I(t)`$ $`=`$ $`{\displaystyle \underset{x=\mathrm{}}{\overset{\mathrm{}}{}}}m(x,t)[1m(x+1,t)]`$ (2.2)
$`S`$ and $`I`$ are representatives of a general class of ”observables” $`M`$ that are sums on $`x`$ of a summand to which each lattice site contributes a factor $`m`$$`1m`$, or $`1`$, i.e., the summand tests for the presence of a specific pattern of visited (”black”) and unvisited (”white”) sites. The following slightly more abstract characterization of the $`M`$ will be needed. Let $`A=\{a\}`$ be a finite set of distinct nonnegative integers $`a`$, such that either $`A=\mathrm{}`$ or, if not, $`A`$ includes the element $`a=0`$. The general observable $`M(t)`$ that we will consider is

$$M(t)=\underset{x=\mathrm{}}{\overset{\mathrm{}}{}}\underset{A}{}\mu _A\underset{aA}{}m(x+a,t)$$
(2.3)
where for $`A=\mathrm{}`$ the product is equal to unity and where $`\{\mu _A\}`$ is a set of numerical coefficients characteristic of $`M`$. When their $`M`$ dependence needs to be indicated we will write $`\mu _A[M]`$. Eqs. (2.2) and (2.1) show that $`S(t)`$ and $`I(t)`$ are of the form of Eq. (2.3) with only two nonzero coefficients, as shown in Table II.

Two further examples of observables of type (2.3) are the total number $`S_1(t)`$ of visited sequences consisting of only a single site, and the total number $`I_1(t)`$ of single-site unvisited sequences. Their coefficients $`\mu _A`$ involve sets $`A`$ of up to three elements; they are easily determined and have also been listed in Table II.

The following remarks, important for later, are verified without much effort. The coefficient $`A_{\mathrm{}}`$ is nonzero if and only if $`M`$ is built up exclusively out of factors $`m`$. Since these correspond to visited sites, that make up the ”bulk” of the support, we will call an $`M`$ of this type a bulk observable. Observables built up exclusively out of factors $`m`$ do not occur, since their expectation value on an infinite lattice is infinite. Hence the remaining observables refer to patterns consisting of both visited and unvisited sites, and we will therefore call them surface observables. \[In the terminology of Refs.  these are ”black” and ”black-and-white” observables. They might also be called ”S-like” and ”I-like”, respectively.\] The distinction between these two subclasses will play a role only in the exponent regime $`1\alpha <2`$.

### 2.2 Basic formulas for averages and correlations

In this work we will first evaluate the $`t\mathrm{}`$ behavior of the averages $`\overline{M(t)}`$. Then we turn to the covariance matrix $`\overline{\mathrm{\Delta }M(t)\mathrm{\Delta }M^{}(t)}`$, where $`M^{}(t)`$ is a second observable with coefficients $`\mu _A^{}`$. Although the authors of Refs.  and  deal with the simple random walk, the larger part of their formal developments also holds for the Riemann walk.

The averages $`\overline{M(t)}`$ and $`\overline{M(t)M^{}(t)}`$ can be obtained as follows . Let $`G(x,t)`$ be the Green function of the one-dimensional Riemann walk, that is, the probability for a walker starting at the origin to occupy site $`x`$ after $`t`$ steps. Let $`\widehat{G}(x,z)=_{t=0}^{\mathrm{}}z^tG(x,t)`$ denote its generating function and let G$`{}_{A}{}^{}(z)`$ be the $`|A|\times |A|`$ matrix of elements $`\widehat{G}(aa^{},z)`$ with $`a,a^{}A`$. From this matrix one constructs the ”inverse sum” $`𝖦_A(z)`$ defined by

$$𝖦_A^1(z)=\underset{a,a^{}A}{}[\text{G}_A^1(z)]_{aa^{}}$$
(2.4)
These scalars satisfy certain elementary relations stated in Appendix A as Properties 1–3. Two functions $`C_M(z)`$ and $`C_{MM^{}}(z)`$ are defined in terms of the $`𝖦_A(z)`$ according to

$`C_M(z)`$ $`=`$ $`{\displaystyle \underset{A\mathrm{}}{}}\mu _A{\displaystyle \frac{1}{𝖦_A(z)}}`$ (2.5)
$`C_{MM^{}}(z)`$ $`=`$ $`{\displaystyle \underset{A\mathrm{}}{}}{\displaystyle \underset{B\mathrm{}}{}}\mu _A\mu _B^{}{\displaystyle \underset{r=\mathrm{}}{\overset{\mathrm{}}{}}}\left[{\displaystyle \frac{1}{𝖦_{A(r+B)}(z)}}{\displaystyle \frac{1}{𝖦_A(z)}}{\displaystyle \frac{1}{𝖦_B(z)}}\right]`$ (2.6)
Here $`A(r+B)`$ denotes the union of the set $`B`$, translated by $`r`$, and $`A`$. The averages $`\overline{M(t)}`$ and $`\overline{M(t)M^{}(t)}`$ are then obtained as

$`\overline{M(t)}`$ $`=`$ $`{\displaystyle \frac{1}{2\pi \text{i}}}{\displaystyle \frac{\text{d}z}{z^{t+1}}\frac{1}{(1z)^2}C_M(z)}`$ (2.7)
$`\overline{M(t)M^{}(t)}`$ $`=`$ $`{\displaystyle \frac{1}{2\pi \text{i}}}{\displaystyle \frac{\text{d}z}{z^{t+1}}\frac{1}{(1z)^2}C_{MM^{}}(z)}`$ (2.8)
where the integrations are counterclockwise around the origin.

The identities (2.4)–(2.8) are fundamental to random walk theory; they hold for any translationally invariant random walk, whether on a finite lattice with periodic boundary conditions or on an infinite lattice. They allow for the calculation, in a very compact way, of many known and new results.

Special cases. When $`M(t)=S(t)`$, the following simplifications occur. The sums on $`A`$ and on $`B`$ in Eqs. (2.5) and (2.6) then have only the single term with $`A=\{0\}`$ and $`B=\{0\}`$, respectively. Furthermore G$`{}_{\{0\}}{}^{}(z)=\widehat{G}(0,z)`$, the matrix G$`{}_{\{0\}(r+\{0\})}{}^{}(z)`$ is two by two, and an easy calculation leads to G$`{}_{\{0\}(r+\{0\})}{}^{}(z)=\frac{1}{2}(\widehat{G}(0,z)+\widehat{G}(r,z))`$. When $`M(t)=I(t)`$, the sum in Eq. (2.5) involves G$`{}_{\{0\}}{}^{}(z)`$ and G$`{}_{\{0,1\}}{}^{}(z)=\frac{1}{2}(\widehat{G}(0,z)+\widehat{G}(1,z))`$. The sums on $`A`$ and $`B`$ in Eq. (2.6) then lead to four terms, which may be evaluated with a little more effort.

### 2.3 Limit $`t\mathrm{}`$ and scaling limit

Explicit evaluation of the general expressions (2.4) -(2.8) is limited in practice by the calculation of the inverse sums $`𝖦_A`$, which require the inversion of a matrix of dimension $`|A|`$. Similarly, evaluation of $`𝖦_{A(r+B)}`$ is an inversion problem of dimension $`|A|+|B|`$ (when $`A`$ and $`r+B`$ have an empty intersection). It turns out that the sum on $`r`$ in Eq. (2.6) can be performed only in the scaling limit

$$z1,|r|\mathrm{}\text{with}\xi =r(1z)^{\frac{1}{\alpha }}\text{fixed}$$
(2.9)
Finally, it will be possible to evaluate the integrals in Eqs. (2.7) and (2.8) only asymptotically for $`t\mathrm{}`$, a limit already implied by Eq. (2.9).

In order to prepare for these limits we rewrite the preceding expressions as follows. Using the simplified notation $`G_0(z)=\widehat{G}(0,z)`$ we split the generating function $`\widehat{G}(x,z)`$ up according to

$$\widehat{G}(x,z)=G_0(z)g(x,z)$$
(2.10)
In full analogy to G$`{}_{A}{}^{}(z)`$ we define g$`{}_{A}{}^{}(z)`$ as the matrix of elements $`g(aa^{},z)`$ with $`a,a^{}A`$, and $`g_A^1`$ as the sum of all elements of g$`{}_{}{}^{1}{}_{A}{}^{}`$. Let now J be the square matrix of elements $`J_{aa^{}}=1`$. Then

$$\text{G}_A(z)=G_0(z)\text{J}\text{g}_A(z)$$
(2.11)
and, by Property 1 of Appendix A,

$$𝖦_A(z)=G_0(z)g_A(z)$$
(2.12)
This splitup will be useful for studying the $`z1`$ behavior of $`𝖦_A(z)`$. Although $`G_0(z)`$ may $`(1\alpha <2)`$ or may not $`(0<\alpha <1)`$ diverge as $`z1`$, the functions $`g(x,z)`$ and $`g_A(z)`$ remain finite in that limit.

We now turn to the inverse sum $`𝖦_{A(r+B)}`$ constructed from the matrix GA∪(r+B). The dimension of this matrix is typically $`|A|+|B|`$. Let JAB be the $`|A|\times |B|`$ matrix with all $`J_{ab}^{AB}=1`$. We then have (for $`AB=\mathrm{}`$)

$$\text{G}_{A(r+B)}(z)=\left(\begin{array}{cc}\text{G}_A(z)& \widehat{G}(r,z)\text{J}^{AB}\\ \widehat{G}(r,z)\text{J}^{BA}& \text{G}_B(z)\end{array}\right)+\left(\begin{array}{cc}0& \text{V}\\ \text{V}^T& 0\end{array}\right)$$
(2.13)
where V is the matrix of elements

$$V_{a,r+b}=\widehat{G}(r+ba,z)\widehat{G}(r,z)aA,bB$$
(2.14)
and VT is its transpose. The first matrix on the RHS of (2.13) has the form (A.4) of Appendix A. Applying Property 3 to that matrix we conclude that

$$\frac{1}{𝖦_{A(r+B)}(z)}=\frac{𝖦_A(z)+𝖦_B(z)2\widehat{G}(r,z)}{𝖦_A(z)𝖦_B(z)\widehat{G}^2(r,z)}+𝒪(\text{V}^2)$$
(2.15)
where we anticipate, and will have to show later, that V is small, that the correction terms are of order $`𝒪(𝐕^2)`$, and that they are negligible for our purpose.

Further analysis depends on the exponent $`\alpha `$. We consider the two main regimes $`0<\alpha <1`$ and $`1<\alpha <2`$ in Secs. 3 and 4, respectively. The exceptional values $`\alpha =\frac{2}{3}`$ and $`\alpha =1`$ are discussed in Sec. 5.

### 2.4 Riemann walk Green function

All quantities of interest have been expressed above in terms of the Riemann walk Green function $`\widehat{G}(x,z)`$. An elementary calculation yields

$`\widehat{G}(x,z)`$ $`=`$ $`{\displaystyle _\pi ^\pi }{\displaystyle \frac{\text{d}q}{2\pi }}{\displaystyle \frac{\text{e}^{\text{i}qx}}{1z\lambda (q)}}`$ (2.16)
$`\lambda (q)`$ $`=`$ $`{\displaystyle \frac{1}{\zeta (1+\alpha )}}{\displaystyle \underset{\mathrm{}=1}{\overset{\mathrm{}}{}}}\mathrm{}^{1\alpha }\mathrm{cos}\mathrm{}q`$ (2.17)
The $`q0`$ behavior of $`\lambda (q)`$ is crucial for the large scale features of the Riemann walk. It is known  that

$$\lambda (q)=1C_\alpha |q|^\alpha +𝒪(q^2)(q0)$$
(2.18)
where for completeness we state the explicit expression

$$C_\alpha ^1=2\zeta (1+\alpha )\mathrm{\Gamma }(1+\alpha )/[\pi \mathrm{sin}(\alpha \pi /2)]$$
(2.19)
One finds by standard methods (see e.g. ) that in the limit $`z1`$ the Green function in the origin $`G_0(z)`$ has the asymptotic expansion

$`G_0(z)`$ $`=`$ $`G_0(1)B_\alpha (1z)^{\frac{1}{\alpha }1}+𝒪(1z)(0<\alpha <1;\alpha \frac{1}{2})`$ (2.20)
$`G_0(z)`$ $`=`$ $`\frac{1}{3}\mathrm{log}[c(1z)^1]+𝒪(1z)(\alpha =1)`$ (2.21)
$`G_0(z)`$ $`=`$ $`A_\alpha (1z)^{1+\frac{1}{\alpha }}+𝒪(1)(1<\alpha <2)`$ (2.22)
where $`B_\alpha `$ and $`A_\alpha `$ are the constants

$`B_\alpha `$ $`=`$ $`C_\alpha ^{1/\alpha }/[2\mathrm{sin}(\pi /\alpha )](\frac{1}{2}<\alpha <1)`$ (2.23)
$`A_\alpha `$ $`=`$ $`1/[2\alpha C_\alpha ^{1/\alpha }\mathrm{sin}(\pi /\alpha )](1<\alpha <2)`$ (2.24)
and $`c`$ is a constant such that there is no $`𝒪(1)`$ term in Eq. (2.21). In Eqs. (2.20)–(2.22) and elsewhere we use the following convention. The symbol $`𝒪(X)`$ indicates terms that are of order X  in the applicable limit ($`X0`$ or $`X\mathrm{}`$); this however is not to say that all preceding terms are larger. Thus, the nonanalytic term in Eq. (2.20) is larger than the $`𝒪(1z)`$ terms only for $`\frac{1}{2}<\alpha <1`$. For $`0<\alpha <\frac{1}{2}`$ it is present only as a correction to the $`𝒪(1z)`$ terms; the expression for its coefficient $`B_\alpha `$ in that regime is different from Eq. (2.23) but will not be needed. For the borderline case $`\alpha =\frac{1}{2}`$, excluded from Eq. (2.20), we have $`G_0(z)=G_0(1)\frac{2}{\pi }(1z)\mathrm{log}(1z)^1+𝒪(1z)`$; but this special nonanalytic behavior will stay subdominant everywhere in the remainder.

From Eqs. (2.16) and (2.18) one deduces that in the scaling limit (2.9)

$$\widehat{G}(r,z)(1z)^{\frac{1}{\alpha }1}F(\xi )(0<\alpha <2)$$
(2.25)
where $`\xi =r(1z)^{\frac{1}{\alpha }}`$ and $`F(\xi )`$ is the scaling function

$$F(\xi )=_{\mathrm{}}^{\mathrm{}}\frac{\text{d}k}{2\pi }\frac{\text{e}^{\text{i}k\xi }}{1+C_\alpha |k|^\alpha }$$
(2.26)
For $`\xi 0`$ it behaves as

$`F(\xi )`$ $``$ $`2\alpha \zeta (1+\alpha )/[\pi \mathrm{sin}(\alpha \pi )]\xi ^{1+\alpha }(0<\alpha <1)`$ (2.27)
$`F(\xi )`$ $``$ $`\frac{1}{3}\mathrm{log}\xi ^1(\alpha =1)`$ (2.28)
$`F(\xi )`$ $`=`$ $`A_\alpha +𝒪(\xi ^{\alpha 1})(1\alpha <2)`$ (2.29)
In Secs. 4 and 5 we will also use the function $`f(r,z)`$ defined by

$$\widehat{G}(r,z)=G_0(z)f(r,z)(1\alpha <2)$$
(2.30)
In the scaling limit one has $`f(r,z)f(\xi )=F(\xi )/F(0)`$ when $`1<\alpha <2`$.

### 2.5 Support at $`t=\mathrm{}`$

Whereas the remainder of this paper deals with the large $`t`$ behavior, we briefly comment here on the structure of the support $`𝒮(t)`$ at $`t=\mathrm{}`$.

As is well-known ), random walks are recurrent (are transient) if $`G_0(1)=\mathrm{}`$ (if $`G_0(1)<\mathrm{}`$). The Riemann walk of this work is recurrent for $`1\alpha <2`$, which means that all sites are visited with probability 1, and that at $`t=\mathrm{}`$ the support $`𝒮(\mathrm{})`$ coincides with the full one-dimensional lattice.

For $`0<\alpha <1`$, however, the Riemann walk is transient, so that at $`t=\mathrm{}`$ the support $`𝒮(\mathrm{})`$ will still be only a subset of the full lattice. We may estimate the average number of visited sites $`\mathrm{\Sigma }_L`$ between $`x=L`$ and $`x=L`$ in $`𝒮(\mathrm{})`$. According to standard random walk theory

$$\mathrm{\Sigma }_L=\underset{x=L}{\overset{L}{}}\frac{\widehat{G}(x,1)}{G_0(1)}$$
(2.31)
Upon substituting (2.16) in (2.31) one easily evaluates $`\mathrm{\Sigma }_L`$ for asymptotically large $`L`$, with the result that $`\mathrm{\Sigma }_LL^\alpha `$. It follows that the support $`𝒮(\mathrm{})`$ has fractal dimension $`d_𝒮=\alpha `$.

## 3 Riemann walk of exponent $`\mathrm{\hspace{0.17em}0}<\alpha <1`$

### 3.1 Averages

The large time behavior of $`\overline{M(t)}`$ comes from the behavior of $`C_M(z)`$, defined in Eq. (2.5), in the limit $`z1`$. From Eq. (2.10) and the explicit expressions (2.16) and (2.17) it may be shown that $`g(x,z)=g(x,1)+𝒪(1z)`$ for all $`0<\alpha <1`$. Hence

$$g_A(z)=g_A(1)+𝒪(1z)$$
(3.1)
after which it follows from Eqs. (2.12), (2.22), and (3.1) that

$$𝖦_A(z)=𝖦_A(1)B_\alpha (1z)^{\frac{1}{\alpha }1}+𝒪(1z)$$
(3.2)
Inverting this relation and substituting in Eq. (2.5) gives

$$C_M(z)=m_1B_\alpha m_2(1z)^{\frac{1}{\alpha }1}B_\alpha ^2m_3(1z)^{\frac{2}{\alpha }2}+\mathrm{}+𝒪(1z)$$
(3.3)
where the $`m_n`$ are determined by the coefficients $`\mu _A`$ of the observable $`M`$ according to

$$m_n[M]=\underset{A\mathrm{}}{}\frac{\mu _A}{𝖦_A^n(1)}(n=1,2,\mathrm{};\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0}<\alpha <1)$$
(3.4)
In Eq. (3.3) the dots stand for a power series in $`(1z)^{\frac{1}{\alpha }1}`$ and the number of nonanalytic terms between the zeroth and the first power of $`1z`$ is equal to $`n_\alpha \frac{\alpha }{1\alpha }1`$. That is, $`n_\alpha `$ is zero for $`0<\alpha <\frac{1}{2}`$ and, as $`\alpha `$ goes up, jumps to $`1,2,3,\mathrm{}`$ at $`\alpha =\frac{1}{2},\frac{2}{3},\frac{3}{4},\mathrm{},`$ respectively. After Laplace inversion we get for the average $`\overline{M(t)}`$ in the limit $`t\mathrm{}`$ the explicit result

$$\overline{M(t)}=m_1t+\frac{B_\alpha }{\mathrm{\Gamma }(3\frac{1}{\alpha })}m_2t^{2\frac{1}{\alpha }}+\frac{B_\alpha ^2}{\mathrm{\Gamma }(4\frac{2}{\alpha })}m_3t^{3\frac{2}{\alpha }}+\mathrm{}+𝒪(1)$$
(3.5)
where the number of nonanalytic terms between the leading and the $`𝒪(1)`$ term is again equal to $`n_\alpha `$.

### 3.2 Correlations

In this subsection we consider two – possibly equal – observables $`M`$ and $`M^{}`$, represented by sets of coefficients $`\{\mu _A\}`$ and $`\{\mu _A^{}\}`$, respectively, and wish to study their correlation. The starting point is Eq. (2.6) for $`C_{MM^{}}(z)`$, in which we substitute Eq. (2.15). Whereas $`𝖦_A(z)`$ and $`𝖦_B(z)`$ tend to finite values in the limit $`z1`$, the Green function $`\widehat{G}(r,z)`$ vanishes in that limit when taken with $`\xi `$ fixed. This suggests that we expand in powers of $`\widehat{G}(r,z)`$,

$$\frac{1}{𝖦_{A(r+B)}(z)}\frac{1}{𝖦_A(z)}\frac{1}{𝖦_B(z)}=\underset{n=1}{\overset{\mathrm{}}{}}C_{AB}^{(n)}(z)\widehat{G}^n(r,z)+𝒪(\text{V}^2)$$
(3.6)
with coefficients

$$C_{AB}^{(n)}(z)=\{\begin{array}{cc}[𝖦_A(z)+𝖦_B(z)][𝖦_A(z)𝖦_B(z)]^{\frac{n}{2}1}\hfill & (n\text{ even})\hfill \\ 2[𝖦_A(z)𝖦_B(z)]^{\frac{n+1}{2}}\hfill & (n\text{ odd})\hfill \end{array}$$
(3.7)
The sum over space that occurs in Eq. (2.6) leads us to now consider the sums $`_r\widehat{G}^n(r,z)`$. From conservation of probability one finds that for $`n=1`$

$$\underset{r}{}\widehat{G}(r,z)=(1z)^1$$
(3.8)
For general $`n`$ the calculation of $`_r\widehat{G}^n(r,z)`$ is slightly more laborious; after substituting Eq. (2.16) for $`\widehat{G}`$ one finds by explicit expansion in powers of $`1z`$ that for $`z1`$ the sum on $`r`$ behaves as

$$\underset{r}{}\widehat{G}^n(r,z)F_{\alpha ,n}(1z)^{1+(n1)(\frac{1}{\alpha }1)}+𝒪(1)(0<\alpha <1;\alpha 1\frac{1}{n})$$
(3.9)
For $`n=1\frac{1}{1\alpha }`$ (with $`n=1,2,\mathrm{}`$) the sum on $`r`$ instead behaves as $`\mathrm{log}(1z)`$; this happens, in particular, for $`n=3`$ when $`\alpha =\frac{2}{3}`$, a case studied separately in Sec. 5.

We will now continue to consider the generic case. The nonanalytic term on the RHS of Eq. (3.9) dominates the $`𝒪(1)`$ term only for $`n<\frac{1}{1\alpha }`$. In that case (3.9) follows just from the scaling form (2.22) of $`\widehat{G}`$ and from the $`\xi 0`$ behavior (2.27) of $`F(\xi )`$. One then finds for the prefactor $`F_{\alpha ,n}`$ the expression

$$F_{\alpha ,n}=2_0^{\mathrm{}}\text{d}\xi F^n(\xi )(1\frac{1}{n}<\alpha <1)$$
(3.10)
For $`0<\alpha <1\frac{1}{n}`$ the expression for $`F_{\alpha ,n}`$ is different and will not be needed. Eq. (3.10) shows that the main contribution to the sum on $`r`$ comes from $`\xi 1`$, that is, from $`r(1z)^{\frac{1}{\alpha }}`$. For $`n=1`$ and $`n=2`$ the integral (3.10) yields the explicit results $`F_{\alpha ,1}=1`$, in agreement with Eq. (3.8), and $`F_{\alpha ,2}=(\frac{1}{\alpha }1)B_\alpha `$, respectively. For $`n>\frac{1}{1\alpha }`$ the sum on $`r`$ draws its main contribution from the short distance (nonscaling) regime $`r1`$, and is of $`𝒪(1)`$ for $`z1`$.

By successively substituting Eq. (3.7) in Eq. (3.6), neglecting the $`𝒪`$(V$`{}_{}{}^{2})`$ terms in that equation – which is justified in Appendix B – , then substituting Eq. (3.6) in Eq. (2.6), expanding $`𝖦_A(z)`$ and $`𝖦_B(z)`$ according to Eq. (3.2), and using the $`z1`$ behavior of $`_r\widehat{G}^n(r,z)`$ obtained in Eqs. (3.8) and (3.9) we find

$`C_{MM^{}}(z)=2(1z)^1m_1m_1^{}`$
$`(1z)^{\frac{1}{\alpha }2}B_\alpha \left(3\frac{1}{\alpha }\right)(m_1m_2^{}+m_2m_1^{})`$
$`(1z)^{\frac{2}{\alpha }3}\left[\left(4\frac{2}{\alpha }\right)B_\alpha ^2(m_1m_3^{}+m_3m_1^{}+m_2m_2^{})+2F_{\alpha ,3}m_2m_2^{}\right]`$
$`\mathrm{}+𝒪\left(1\right)`$ (3.11)
Here the dots stand for terms of order $`(1z)^{1+k(\frac{1}{\alpha }1)}`$, with $`k=3,4,\mathrm{};`$ and the $`m_n^{}`$ are related to $`M^{}`$ in the same way as the $`m_n`$ are to $`M`$. For $`t\mathrm{}`$ we therefore find by substituting (3.11) in Eq. (2.8) and Laplace inverting

$$\begin{array}{cc}& \overline{M(t)M^{}(t)}=m_1m_1^{}t^2\hfill \\ & +\frac{B_\alpha }{\mathrm{\Gamma }(3\frac{1}{\alpha })}(m_1m_2^{}+m_2m_1^{})t^{3\frac{1}{\alpha }}\hfill \\ & +\left[\frac{B_\alpha ^2}{\mathrm{\Gamma }(4\frac{2}{\alpha })}(m_1m_3^{}+m_3m_1^{}+m_2m_2^{})+\frac{2F_{\alpha ,3}}{\mathrm{\Gamma }(5\frac{2}{\alpha })}m_2m_2^{}\right]t^{4\frac{2}{\alpha }}\hfill \\ & +\mathrm{}+𝒪(t)\hfill \end{array}$$

(3.12)
The successive terms in the above series all have one power of $`t`$ more than the corresponding terms in the series (3.5) for $`\overline{M(t)}`$, and the number of nonanalytic terms between the leading and the $`𝒪(t)`$ term is once more equal to $`n_\alpha 1`$. The product $`\overline{M(t)}`$ $`\overline{M^{}(t)}`$, which follows from Eq. (3.5), now has to be subtracted from the series (3.2). This exactly cancels the terms in (3.2) proportional to $`t^2`$ and to $`t^{3\frac{1}{\alpha }}`$ but leaves those proportional to $`t^{4\frac{2}{\alpha }}`$ and of $`𝒪(t)`$. The $`t^{4\frac{2}{\alpha }}`$ terms are leading only if $`\frac{2}{3}<\alpha <1`$. Hence we find for the correlation between observables $`M`$ and $`M^{}`$ in the limit $`t\mathrm{}`$

$$\overline{\mathrm{\Delta }M(t)\mathrm{\Delta }M^{}(t)}_\alpha ^2m_2m_2^{}t^{4\frac{2}{\alpha }}$$
(3.13)
valid for $`\frac{2}{3}<\alpha <1`$, and in which

$$_\alpha ^2=\frac{2F_{\alpha ,3}}{\mathrm{\Gamma }(5\frac{2}{\alpha })}+\frac{B_\alpha ^2}{\mathrm{\Gamma }(4\frac{2}{\alpha })}\frac{B_\alpha ^2}{\mathrm{\Gamma }^2(3\frac{1}{\alpha })}$$
(3.14)
We have supposed here that $`m_2,m_2^{}0`$. The preceding analysis changes when either of these two coefficients vanishes. We do not know of any physically interesting examples where this happens, and do not pursue our analysis in this direction.

The borderline case $`\alpha =\frac{2}{3}`$ is considered in Sec. 5. In the interval $`0<\alpha <\frac{2}{3}`$ the calculation of the present section applies, but with the result that

$$\overline{\mathrm{\Delta }M(t)\mathrm{\Delta }M^{}(t)}\kappa _{MM^{}}t(0<\alpha <\frac{2}{3})$$
(3.15)
in which the coefficient $`\kappa _{MM^{}}`$ has contributions from the $`𝒪(1)`$ terms in Eq. (3.5) and the $`𝒪(t)`$ terms in Eq. (3.2), and does not factor into an $`M`$ and an $`M^{}`$ dependent constant. This difference between Eqs. (3.13) and (3.15) is crucial for the phenomenon of universality discussed in Sec. 6.

## 4 Riemann walk of exponent $`\mathrm{\hspace{0.17em}1}<\alpha <2`$

### 4.1 Averages

The calculation of $`\overline{M(t)}`$ starts again from the series (2.5) for $`C_M(z)`$. The calculation in the exponent regime $`1<\alpha <2`$ is different from that of the preceding section because now $`G_0(z)`$ diverges as $`z1`$. Since $`g_A(z)`$ remains finite for $`z1`$, this suggests that we use Eq. (2.12) and expand $`𝖦_A(z)`$ in powers of $`g_A(z)/G_0(z)`$. This yields

$$C_M(z)=\underset{A\mathrm{}}{}\mu _A\frac{1}{G_0(z)}\left[\mathrm{\hspace{0.17em}\hspace{0.17em}1}+\frac{g_A(z)}{G_0(z)}+\frac{g_A^2(z)}{G_0^2(z)}+𝒪(\frac{g_A^3}{G_0^3})\right]$$
(4.1)
We now substitute in Eq. (4.1) the expansion (2.22) for $`G_0(z)`$ and use that $`g_A(1)`$ is finite. The result is a power series in $`1z`$ in which there appear coefficients that we denote again by $`m_n`$ but that are defined for $`1\alpha <2`$ as

$$m_n[M]=\underset{A\mathrm{}}{}\mu _Ag_A^n(1)(n=0,1,2,\mathrm{};\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}1}\alpha <2)$$
(4.2)
It will turn out that we need only the leading term, which is

$$C_M(z)\{\begin{array}{c}A_\alpha ^1m_0(1z)^{1\frac{1}{\alpha }}(m_00)\hfill \\ A_\alpha ^2m_1(1z)^{2\frac{2}{\alpha }}(m_0=0,m_10)\hfill \end{array}$$
(4.3)
We pause to note that in the terminology of Sec. 2.1 the condition $`m_00`$ characterizes the bulk or ”S-like” observables, and the condition $`m_0=0`$ the surface or ”I-like” observables. This is the first equation where a difference appears between these two subclasses; in its analog, Eq. (3.5) of the preceding section, no such distinction appears.

Upon using Eq. (4.3) in Eq. (2.7) we obtain after an inverse Laplace transformation the asymptotic expansion of $`\overline{M(t)}`$ as $`t\mathrm{}`$,

$$\overline{M(t)}\{\begin{array}{c}[A_\alpha \mathrm{\Gamma }(1+\frac{1}{\alpha })]^1m_0t^{\frac{1}{\alpha }}(m_00)\hfill \\ [A_\alpha \mathrm{\Gamma }(\frac{2}{\alpha })]^1m_1t^{\frac{2}{\alpha }1}(m_0=0,m_10)\hfill \end{array}$$
(4.4)
where the dots indicate terms of lower order in $`t`$.

### 4.2 Correlations

For the calculation of the correlation $`\overline{M(t)M^{}(t)}`$ via Eqs. (2.8) and (2.6) we have to return again to expression (2.15) for $`1/𝖦_{A(r+B)}(z)`$, which is needed in Eq. (2.6). We use Eqs. (2.15), (2.12), and (2.30) to rewrite this quantity as

$$\frac{1}{𝖦_{A(r+B)}(z)}=\frac{1}{G_0(z)}\frac{2(1f(r,z))\frac{g_A(z)}{G_0(z)}\frac{g_B(z)}{G_0(z)}}{1f^2(r,z)\frac{g_A(z)}{G_0(z)}\frac{g_B(z)}{G_0(z)}+\frac{g_A(z)g_B(z)}{G_0^2(z)}}+𝒪(\text{V}^2)$$
(4.5)
The function $`f(r,z)`$ was defined in Eq. (2.30). In the scaling limit $`f(r,z),`$ $`g_A(z)`$, and $`g_B(z)`$ have finite limits, whereas $`G_0(z)`$ diverges. An expansion in inverse powers of $`G_0(z)`$ corresponds therefore to an expansion in ascending powers of $`1z`$. Writing for short $`f,g_A,g_B`$, and $`G_0`$ when $`f(r,z),`$$`g_A(z),`$$`g_B(z)`$, and $`G_0(z)`$ are meant, we find after a straightforward calculation

$`C_{MM^{}}(z)`$ $`=`$ $`{\displaystyle \underset{r=\mathrm{}}{\overset{\mathrm{}}{}}}{\displaystyle \frac{1}{G_0}}{\displaystyle \frac{f}{1+f}}{\displaystyle \underset{A,B\mathrm{}}{}}\mu _A\mu _B^{}[\mathrm{\hspace{0.17em}\hspace{0.17em}2}+(2+f){\displaystyle \frac{g_A+g_B}{G_0}}`$ (4.6)
$`+{\displaystyle \frac{22ff^3}{(1f)(1+f)^2}}{\displaystyle \frac{g_A^2+g_B^2}{G_0^2}}`$
$`+{\displaystyle \frac{2}{(1f)(1+f)^2}}{\displaystyle \frac{g_Ag_B}{G_0^2}}+𝒪({\displaystyle \frac{g_A^3}{G_0^3}},{\displaystyle \frac{g_B^3}{G_0^3}})]`$
Let us write $`m_n^{}=m_n[M^{}]`$ for the coefficients that characterize the observable $`M^{}`$. The two distinct cases described by Eq. (4.4) now lead to the following possibilities.

Case (i): $`m_00`$ and $`m_0^{}0`$. In this case the leading term in the expression in brackets in Eq. (4.6) survives under the sum on $`A`$ and $`B`$.

Case (ii): $`m_0=0,m_10`$, and $`m_0^{}0`$. In this case in order to survive a term in the bracketed expression should contain at least one factor $`g_A(z)`$.

Case (iii): $`m_0=m_0^{}=0`$ but $`m_10`$ and $`m_1^{}0`$. In this case a term in order to survive must contain at least one factor $`g_A(z)`$ and one factor $`g_B(z)`$.
Upon using in each of these cases for $`G_0(z)`$the expansion (2.22), passing to the scaling limit, and writing $`f`$ for $`f(\xi )`$, we find that the result is

$$C_{MM^{}}(z)\{\begin{array}{c}A_\alpha ^1(1z)^{1\frac{2}{\alpha }}f_{00}m_0m_0^{}\hfill \\ A_\alpha ^1(1z)^{2\frac{3}{\alpha }}f_{10}m_1m_0^{}\hfill \\ A_\alpha ^1(1z)^{3\frac{4}{\alpha }}f_{11}m_1m_1^{}\hfill \end{array}$$
(4.7)
in the three cases (i), (ii), and (iii), respectively; here the coefficients $`f_k\mathrm{}`$ represent the integrals

$`f_{00}`$ $`=`$ $`4{\displaystyle _0^{\mathrm{}}}\text{d}\xi f(1+f)^1`$
$`f_{10}`$ $`=`$ $`2{\displaystyle _0^{\mathrm{}}}\text{d}\xi f(2+f)(1+f)^1`$ (4.8)
$`f_{11}`$ $`=`$ $`4{\displaystyle _0^{\mathrm{}}}\text{d}\xi f(1f)^1(1+f)^3`$
After substituting Eqs. (4.7) in Eq. (2.8) and carrying out the inverse Laplace transformation we find, in the limit $`t\mathrm{}`$,

$$\overline{M(t)M^{}(t)}\{\begin{array}{cc}\hfill \mathrm{\Gamma }^1(1+\frac{2}{\alpha })& A_\alpha ^1f_{00}m_0m_0^{}t^{\frac{2}{\alpha }}\hfill \\ \hfill \mathrm{\Gamma }^1(\frac{3}{\alpha })& A_\alpha ^2f_{10}m_1m_0^{}t^{\frac{3}{\alpha }1}\hfill \\ \hfill \mathrm{\Gamma }^1(1+\frac{4}{\alpha })& A_\alpha ^3f_{11}m_1m_1^{}t^{\frac{4}{\alpha }2}\hfill \end{array}$$
(4.9)
respectively, for the three cases distinguished above. Upon combining these results with those of Section 4.1 one obtains, for $`t\mathrm{}`$,

$$\overline{\mathrm{\Delta }M(t)\mathrm{\Delta }M^{}(t)}\{\begin{array}{c}_{\alpha 00}^2m_0m_0^{}t^{\frac{2}{\alpha }}\hfill \\ _{\alpha 10}^2m_1m_0^{}t^{\frac{3}{\alpha }1}\hfill \\ _{\alpha 11}^2m_1m_1^{}t^{\frac{4}{\alpha }2}\hfill \end{array}$$
(4.10)
in which the coefficients $`_{\alpha k\mathrm{}}`$ are given by

$$\begin{array}{c}_{\alpha 00}^2=A_\alpha ^1[f_{00}\mathrm{\Gamma }^1(1+\frac{2}{\alpha })A_\alpha ^1\mathrm{\Gamma }^2(1+\frac{1}{\alpha })]\hfill \\ _{\alpha 10}^2=A_\alpha ^2[f_{10}\mathrm{\Gamma }^1(\frac{3}{\alpha })A_\alpha ^1\mathrm{\Gamma }^1(1+\frac{1}{\alpha })\mathrm{\Gamma }^1(\frac{2}{\alpha })]\hfill \\ _{\alpha 11}^2=A_\alpha ^3[f_{11}\mathrm{\Gamma }^1(1+\frac{4}{\alpha })A_\alpha ^1\mathrm{\Gamma }^2(\frac{2}{\alpha })]\hfill \end{array}$$
(4.11)
in the three cases (i), (ii), and (iii) defined above, respectively. We recall that the $`f_k\mathrm{}`$ on the RHS of Eq. (4.11) are given by Eq. (4.8) as integrals on $`f(\xi )`$, with $`f(\xi )`$ in turn given by Eq. (2.26).

## 5 Riemann walk of exponents $`\alpha =\frac{2}{3}`$ and $`\alpha =1`$

### 5.1 Exponent $`\alpha =\frac{2}{3}`$

In this special case $`\overline{M(t)}`$ is still given by Eq. (3.5). However, the calculation of $`\overline{M(t)M^{}(t)}`$ has to be reconsidered, as signalled by the fact that $`F_{\alpha ,3}`$ in Eq. (3.11) diverges for $`\alpha \frac{2}{3}^+`$. In order to calculate $`_r\widehat{G}(r,z)`$ we cannot now use Eq. (3.9). Instead we replace $`\widehat{G}(r,z)`$ by its scaling form (2.22) but take into account that $`|\xi |`$ has a lower cutoff $`|\xi |\text{cst}\times (1z)^{\frac{2}{3}}`$. This gives

$`{\displaystyle \underset{r}{}}\widehat{G}^3(r,z)`$ $``$ $`2{\displaystyle _{\text{cst}\times (1z)^{2/3}}^{\mathrm{}}}\text{d}\xi F^3(\xi )`$ (5.1)
$``$ $`\frac{1}{2}𝒞^2\mathrm{log}(1z)^1+𝒪(1)(z1)`$
where in the second step we used Eq. (2.24) and found for the coefficient the value $`𝒞^2=2^{12}3^{11/2}\pi ^3\zeta ^3(\frac{3}{2})`$. In this case, due to the $`\mathrm{log}(1z)`$ in the equation above, $`\overline{M(t)M^{}(t)}`$ is larger than the product $`\overline{M(t)}`$ $`\overline{M^{}(t)}`$ by a factor $`\mathrm{log}t`$, and determines by itself alone the final result, which reads

$$\overline{\mathrm{\Delta }M(t)\mathrm{\Delta }M^{}(t)}𝒞^2m_2m_2^{}t\mathrm{log}t(t\mathrm{})$$
(5.2)
This $`t\mathrm{log}t`$ behavior is the same as in the well-known case of the simple random walk in spatial dimension $`d=3`$ .

### 5.2 Exponent $`\alpha =1`$

The case of Lévy exponent $`\alpha =1`$ is subtler than the others. Since it is closely analogous to the simple random walk in dimension $`d=2`$ , we will not present all steps in detail. Eq. (2.21) shows that for $`\alpha =1`$ the Green function in the origin, $`G_0(z)`$, diverges as $`z1`$. We can therefore expand $`C_M(z)`$ as a series in the same way as in Eq. (4.1). Since here again $`g_A(z)=g_A(1)+𝒪(1z)`$, and in view of the logarithmic behavior (2.21), this series now leads to an expansion of $`C_M(z)`$ in inverse powers of $`G_0(z)`$. If the first nonzero term is of order $`k+1`$, then we have explicitly

$$C_M(z)=m_kG_0^{k1}(z)m_{k+1}G_0^{k2}(z)m_{k+2}G_0^{k3}(z)+\mathrm{}$$
(5.3)
with the $`m_n`$ defined by Eq. (4.2). The cases of physical interest have $`k=0`$ (bulk observables) or $`k=1`$ (surface observables), but it will be notationally convenient to keep $`k`$ as a parameter. We will also refer to it as the order of $`M`$.

To find $`C_{MM^{}}(z)`$ we may still start from Eq. (4.5), but now the expansion of this equation runs differently. The reason is that for $`z1`$ at fixed $`\xi `$ the function $`f(r,z)`$ (defined by (2.30)) behaves as $`F(\xi )/G_0(z)`$ and so is of the same order as $`g_A(z)/G_0(z)`$ and $`g_B(z)/G_0(z)`$. We therefore have to perform a double expansion of the RHS of Eq. (4.5) in terms of on the one hand $`F/G_0`$ and on the other hand $`g_A/G_0`$ and $`g_B/G_0`$. The sum on $`r`$, which in the scaling limit becomes an integral on $`\xi `$, then leads to the appearance of coefficients $`F_{1,n}`$ defined as in Eq. (3.10) but with $`\alpha =1`$. Special cases are $`F_{1,1}=1`$ and $`F_{1,2}=\frac{1}{3}`$. Let $`m_k`$ and $`m_k^{}^{}`$ be the first nonzero coefficients in the expansions of $`C_M(z)`$ and $`C_M^{}(z)`$, respectively. Then we find for $`C_{MM^{}}(z)`$, retaining only the three leading order terms in the limit $`z1`$,

$`C_{MM^{}}(z)`$ $`{\displaystyle \frac{1}{(1z)G_0^{k+k^{}+2}(z)}}[2m_km_k^{}^{}`$
$``$ $`G_0^1(z)\left(\frac{1}{3}(k+k^{}+2)a_2m_km_k^{}^{}2(m_km_{k^{}+1}^{}+m_{k+1}m_k^{}^{})\right)`$
$`+`$ $`G_0^2(z)(2F_{1,3}(k+1)(k^{}+1)a_3m_km_k^{}^{}`$
$`\frac{1}{3}(k+k^{}+3)a_2(m_km_{k^{}+1}^{}+m_{k+1}m_k^{}^{})`$
$`+2(m_km_{k^{}+2}^{}+m_{k+1}m_{k^{}+1}^{}+m_{k+2}m_k^{}^{}))]`$ (5.4)
The inverse Laplace transforms of $`C_M(z)`$ and $`C_{MM^{}}(z)`$ may be found with the help of the explicit expression (2.21) for $`G_0(z)`$ and the integrals of Ref. . We state only the explicit result for $`\overline{M(t)}`$, which is, for $`t\mathrm{}`$,

$`\overline{M(t)}{\displaystyle \frac{3^{k+1}t}{\mathrm{log}^{k+1}ct}}[`$ $`m_k+{\displaystyle \frac{1}{\mathrm{log}ct}}\left((1\gamma )(k+1)m_k+3m_{k+1}\right)`$
$`+{\displaystyle \frac{1}{\mathrm{log}^2ct}}((1\frac{1}{12}\pi ^2\gamma +\frac{1}{2}\gamma ^2)(k+1)(k+2)m_k`$
$`3(1\gamma )(k+2)m_{k+1}+9m_{k+2})]`$ (5.5)
in which $`\gamma =0.577215\mathrm{}`$ denotes Euler’s constant. Both $`\overline{M(t)M^{}(t)}`$ and the product $`\overline{M(t)}`$ $`\overline{M^{}(t)}`$ the appear as $`t^2`$ times a power series in $`1/\mathrm{log}ct`$ of which the leading term is of order $`k+k^{}+2`$, and in which the three leading orders have to be retained. Upon carrying out the subtraction one finds that the two leading orders cancel and the correlation $`\overline{\mathrm{\Delta }M(t)\mathrm{\Delta }M^{}(t)}`$ appears to be proportional to $`t^2/\mathrm{log}^{k+k^{}+4}ct`$. Explicitly, as $`t\mathrm{}`$,

$$\overline{\mathrm{\Delta }M(t)\mathrm{\Delta }M^{}(t)}𝒜^2(k+1)(k^{}+1)m_km_k^{}^{}\frac{3^{k+k^{}+2}t^2}{\mathrm{log}^{k+k^{}+4}ct}(t\mathrm{})$$
(5.6)
in which the coefficient $`𝒜`$ is given by

$$𝒜^2=1+(F_{1,3}\frac{1}{6})\pi ^2$$
(5.7)
Numerical evaluation gives $`F_{1,3}=0.27415\mathrm{}`$, whence $`𝒜^2=2.0608\mathrm{}`$. Eq. (5.6) is the same as for the two-dimensional simple random walk , but with a different constant $`𝒜`$.

## 6 Universality of fluctuations

We consider in this section the normalized deviations from average

$$\theta _M(t)=\frac{\mathrm{\Delta }M(t)}{\overline{\mathrm{\Delta }M^2(t)}^{1/2}}$$
(6.1)
These random functions of time satisfy by construction

$$\overline{\theta _M(t)}=0,\overline{\theta _M^2(t)}=1$$
(6.2)
We consider now two arbitrary observables $`M`$ and $`M^{}`$. When $`\frac{2}{3}\alpha <1`$ we have from Eq. (6.1) together with either Eq. (5.2) or Eq. (3.13) that

$$\overline{\theta _M(t)\theta _M^{}(t)}=1$$
(6.3)
It then follows from Eqs. (6.2) and (6.3) that the difference $`\theta _M\theta _M^{}`$ is a random variable of zero average and zero variance. Such a random variable can only be itself equal to zero. We therefore deduce that, when $`\frac{2}{3}\alpha <1`$, in the limit $`t\mathrm{}`$ all $`\theta _M(t)`$ are equal to a single random variable, which we will call $`\mathrm{\Theta }_\alpha (t)`$, thus indicating explicitly its $`\alpha `$ dependence.

When $`1\alpha <2`$ we have for two observables $`M`$ and $`M^{}`$ whose orders, $`k`$ and $`k^{}`$, are equal from Eq. (6.1) and either Eq. (5.6) or Eq. (4.10) again the result (6.3). Hence, when $`1\alpha <2`$, in the limit $`t\mathrm{}`$ all $`\theta _M(t)`$ with $`k=0`$ are equal to a single random variable – that we will call $`\mathrm{\Theta }_{\alpha 0}(t)`$ –, and all $`\theta _M(t)`$ with $`k=1`$ are similarly equal to a single random variable – that we will call $`\mathrm{\Theta }_{\alpha 1}(t)`$. The variables $`\mathrm{\Theta }_\alpha (t)`$, $`\mathrm{\Theta }_{\alpha 0}(t)`$, and $`\mathrm{\Theta }_{\alpha 1}(t)`$ are universal in the sense that they are independent of the observables $`M`$ (but depend at most on their order).

In each of these cases the key ingredient necessary for arriving at Eq. (6.3) is the factorization of $`\overline{M(t)M^{}(t)}`$ into an $`M`$ and an $`M^{}`$ dependent part. This also explains why for $`\alpha <\frac{2}{3}`$ the same reasoning fails.

The cross correlation between $`\mathrm{\Theta }_{\alpha 0}(t)`$ and $`\mathrm{\Theta }_{\alpha 1}(t)`$ is easily found from the correlation between a $`\theta _M(t)`$ and a $`\theta _M^{}(t)`$ with $`k=0`$ and $`k^{}=1`$, and use of the second one of Eqs. (4.10). The answer is independent of the choice of $`M`$ and $`M^{}`$, as it had to be, and reads

$$\overline{\mathrm{\Theta }_{\alpha 0}(t)\mathrm{\Theta }_{\alpha 1}(t)}=_{\alpha 10}^2(_{\alpha 00}_{\alpha 11})^1$$
(6.4)
The coefficient ratio on the RHS of this equation depends only on the exponent $`\alpha `$ and must necessarily be less than unity. In the limit $`\alpha 1^+`$ it approaches unity and for $`\alpha <1`$ the distinction between surface and bulk observables is no longer reflected in the fluctuations.

Upon combining all these conclusions we get explicitly

$$\mathrm{\Delta }M(t)=\{\begin{array}{cc}m_2𝒞(t\mathrm{log}t)^{\frac{1}{2}}\mathrm{\Theta }_{\frac{2}{3}}(t)\hfill & (\alpha =\frac{2}{3})\hfill \\ m_2_\alpha t^{2\frac{1}{\alpha }}\mathrm{\Theta }_\alpha (t)\hfill & (\frac{2}{3}<\alpha <1)\hfill \\ m_k𝒜(k+1)t(\frac{1}{3}\mathrm{log}ct)^{k2}\mathrm{\Theta }_1(t)\hfill & (\alpha =1;k=0,1)\hfill \\ m_0_{\alpha 00}t^{\frac{1}{\alpha }}\mathrm{\Theta }_{\alpha 0}(t)\hfill & (1<\alpha <2;k=0)\hfill \\ m_1_{\alpha 11}t^{\frac{2}{\alpha }1}\mathrm{\Theta }_{\alpha 1}(t)\hfill & (1<\alpha <2;k=1)\hfill \end{array}$$
(6.5)
Here all $`M`$ dependence is contained in the coefficients $`m_n`$.

## 7 Conclusions

We have studied a large class of properties $`M(t)`$ of the support of the one-dimensional $`t`$ step Riemann walk. These include the number $`S(t)`$ of distinct sites visited, and the number $`I(t)`$ of sequences of visited sites. The $`M(t)`$ fall into two classes, the bulk or S-like properties, and the surface or I-like properties. The asymptotic laws found in the preceding sections for the averages, variances, and correlation of $`S(t)`$ and $`I(t)`$ have been summarized in Table I in the Introduction.

It appears from that table that in the exponent regime $`0<\alpha 1`$ the ratios $`\overline{\mathrm{\Delta }S^2(t)}^{1/2}/\overline{S(t)}`$ and $`\overline{\mathrm{\Delta }I^2(t)}^{1/2}/\overline{I(t)}`$ tend to zero when $`t\mathrm{}`$, which indicates that the distributions of $`S(t)`$ and of $`I(t)`$ become infinitely narrowly peaked around their average. Hence in this exponent regime the ratio $`s(t)\overline{S(t)}/\overline{I(t)}`$ represents the average number of sites per visited sequence. When $`\alpha `$ is strictly less than unity we have explicitly

$$\underset{t\mathrm{}}{lim}s(t)=\frac{m_1[S]}{m_1[I]}=\frac{\widehat{G}(0,1)\widehat{G}(1,1)}{\widehat{G}(0,1)+\widehat{G}(1,1)}(0<\alpha <1)$$
(7.1)
The first equality is based on Eq. (3.5) and in the second one we used the definition (3.4) of the $`m_n[M]`$ and the remarks at the end of Sec. 2.2. The finiteness of the result (7.1) means that in the large $`t`$ limit every new step of the walk creates a new visited sequence with a finite nonzero probability. This explains that in this exponent regime the asymptotic power laws do not distinguish between bulk and surface properties. For $`\alpha 1^{}`$ expression (7.1) diverges, and when $`\alpha =1`$ the ratio $`s(t)`$ increases logarithmically with $`t`$.

In the exponent regime $`1<\alpha <2`$ the appropriately scaled distributions of $`S(t)`$ and $`I(t)`$ are of finite width even in the limit $`t\mathrm{}`$. The support has an ”interior”, bulk and surface properties have different asymptotic power laws, and $`s(t)\mathrm{}`$ as $`t\mathrm{}`$. In the terminology of critical phenomena, this regime is fluctuation dominated. In this regime the universality of fluctuations holds in a slightly weaker but at least as interesting a sense as for $`0<\alpha 1`$. To describe the fluctuations, not a single but two universal stochastic variables are needed, one applying to the bulk and the other to the surface properties. These two variables become fully correlated in the limit $`\alpha 1^+`$.

Finally we remark that when $`2/\alpha `$ is equal to one of the integers $`2,3,4,\mathrm{}`$, the asymptotic laws of Table I coincide with the ones known to hold for the simple random walk on a lattice of dimension $`d=2/\alpha `$. Similarly, the universality properties for those $`\alpha `$ values have their analogs in the $`d`$-dimensional simple random walk. Hence the ”rule of the effective dimensionality”, which states the correspondence $`\alpha 2/d`$, applies to all properties that we have studied. Of course it must break down when the comparison between the Riemann walk and the simple random walk is refined sufficiently. Also, we have not considered the borderline value $`\alpha =2`$, which is special , and for which this rule fails.

## Acknowledgments

The authors acknowledge support from the French-Brazilian scientific cooperation project CAPES/COFECUB 229/97.

## Appendix A Relations for the inverse sums $`𝖦_A`$ and $`g_A`$

We collect here some elementary matrix algebra relations useful for dealing with the inverse sums $`𝖦_A`$ and $`g_A`$ occurring in the main text. The $`z`$ dependence of these quantities plays no role. The presentation and notation are independent of the body of the paper.

Let $`L`$ be an invertible $`\mathrm{}\times \mathrm{}`$ matrix. We define the ”inverse sum” $`(L)`$ by

$$^1(L)=\underset{i,j}{}L_{ij}^1$$
(A.1)
In the remainder $`\alpha ,\beta ,`$ and $`\gamma `$ will denote constants.

Property 1. Let $`J`$ be the $`\mathrm{}\times \mathrm{}`$ matrix with all $`J_{ij}=1`$, and let $`M`$ be an invertible $`\mathrm{}\times \mathrm{}`$ matrix. Let $`L=\alpha J+\gamma M`$. Then

$$(L)=\alpha +\gamma (M)$$
(A.2)
The proof of this relation is given in Ref. .

Property 2. Let $`M`$ and $`N`$ be invertible matrices of dimensions $`m\times m`$ and $`n\times n`$, respectively, and let $`L`$ be the block diagonal $`\mathrm{}\times \mathrm{}`$ matrix with blocks $`M`$ and $`N`$. Then

$$\frac{1}{(L)}=\frac{1}{(M)}+\frac{1}{(N)}$$
(A.3)
This follows directly from the definition (A.1). The calculation of $`(L)`$ for an $`\mathrm{}\times \mathrm{}`$ matrix may be reduced to an inversion problem of dimension less than $`\mathrm{}`$ also in certain cases where $`L`$ is not block diagonal, as shown below.

Property 3. Let $`J^{mn}`$ be the $`m\times n`$ matrix with all elements equal to 1. Let $`L`$ be $`\mathrm{}\times \mathrm{}`$ and of the form

$$L=\left(\begin{array}{cc}\gamma M& \beta J^{mn}\\ \beta J^{nm}& \gamma N\end{array}\right)$$
(A.4)
Then

$$(L)=\frac{\gamma ^2(M)(N)\beta ^2}{\gamma (M)+\gamma (N)2\beta }$$
(A.5)
To prove this we rewrite $`L`$ as $`L=\beta J+\stackrel{~}{L}`$, where $`J`$ is as before and where

$$\stackrel{~}{L}=\left(\begin{array}{cc}\gamma M\beta J^{mm}& 0\\ 0& \gamma N\beta J^{nn}\end{array}\right)$$
(A.6)
From Property 1 we have that $`(L)=\beta +(\stackrel{~}{L})`$, after which by applying Property 2 and once more Property 1, we obtain after some rearrangement Eq. (A.5). For $`\beta =0`$ Eq. (A.5) reduces to Property 2.

In this work the need for Properties 1 and 3 arises when the limit $`\gamma 0`$ has to be taken. For $`\gamma =0`$ the matrices $`L`$ that occur on the LHS of Eqs. (A.2) and (A.5) are no longer invertible, but these properties allow nevertheless $`(L)`$ to be calculated in that limit.

## Appendix B Corrections to scaling

In Eq. (3.6) we have neglected the $`𝒪`$(V2) terms that appear in Eq. (2.15). Since in the last step that led to Eq. (3.13) the leading order in $`1z`$ went down due to cancellations, we must now check that the $`𝒪`$(V2) terms remain subdominant. In this Appendix we will write Eq. (2.13) in the simplified notation $`\text{G}_{A(r+B)}=\text{G}+\text{W}`$ where G and W are the first and second matrix, respectively, on the RHS of Eq. (2.13). Upon writing the inverse G$`{}_{}{}^{1}{}_{A(r+B)}{}^{}`$ as a perturbation series in W and applying Eq. (2.4) one finds

$$𝖦_{A(r+B)}^1=\underset{\mathrm{}=0}{\overset{\mathrm{}}{}}(1)^{\mathrm{}}\underset{c,c^{}}{}[\text{G}^1(\text{WG}^1)^{\mathrm{}}]_{cc^{}}$$
(B.1)
The $`\mathrm{}=0`$ term of this series is the term shown explicitly on the RHS of Eq. (2.15), and has been the object of study in Sec. 3.2. We will show here that the terms with $`\mathrm{}1`$ produce, in the scaling limit, only higher order corrections to the final result. To this end we first consider $`C_{MM^{}}(z)`$ defined by Eq. (2.6). Let $`R_{\mathrm{}}(z)`$ denote the contribution to $`C_{MM^{}}(z)`$ from the $`\mathrm{}`$th term in Eq. (B.1). In order to estimate the order in $`\mathrm{\hspace{0.17em}1}z`$ of $`R_{\mathrm{}}(z)`$ as $`z1`$ we first deduce from Eqs. (2.14) and (2.25) that in the scaling limit the matrix elements of W behave as $`V_{a,r+b}(1z)^{\frac{2}{\alpha }1}(ba)F^{}(\xi )`$, and that summing on $`r`$ amounts to applying $`(1z)^{\frac{1}{\alpha }}\text{d}\xi `$. This yields the asymptotic proportionality

$$R_{\mathrm{}}(z)(1z)^{\frac{1}{\alpha }}\text{d}\xi \frac{1}{G_0(z)}\left[\frac{(1z)^{\frac{2}{\alpha }1}F^{}(\xi )}{G_0(z)}\right]^{\mathrm{}}$$
(B.2)
where $`G_0(z)`$ represents the order in $`\mathrm{\hspace{0.17em}1}z`$ of the matrix G. When $`\mathrm{}`$ is odd, this integral vanishes by symmetry, which shows that the leading correction is of order V2, as anticipated. For $`\xi 0`$ we have, in virtue of Eq. (2.27), that $`F^{}(\xi )\xi ^{2+\alpha }`$. Hence the $`\xi `$ integral in Eq. (B.2) diverges in the origin for all $`\mathrm{}1`$ when $`0<\alpha <1`$. This signals that the main contribution comes from $`r`$ values near the origin. The order in $`\mathrm{\hspace{0.17em}1}z`$ of $`R_{\mathrm{}}(z)`$ may then be estimated by introducing in the integral the cutoff $`|\xi |\text{cst}\times (1z)^{\frac{1}{\alpha }}`$, which leads to $`R_{\mathrm{}}(z)(1z)^0`$. When $`\mathrm{}1`$ these additive corrections to $`C_{MM^{}}(z)`$ in Eq. (3.11) are negligible, therefore, with respect to the $`(1z)^{\frac{2}{\alpha }3}`$ term which, in the relevant exponent regime $`\frac{2}{3}<\alpha <1`$, determines the final result.


format: text

text:
# 1 Introduction

## 1 Introduction

Lattice QCD simulations are usually performed with periodic boundary conditions (BC). However, other type of BC may be important in certain cases. For example, a comparison between systems with different BC can be used to understand finite size-effects in lattice QCD. Some years ago? $`C`$-periodic BC were studied as an alternative to periodic conditions. Then they were considered with the general idea of studying the spontaneous symmetry breaking aspects of the QCD dynamics in a simple way.? In that work an analysis was done in the continuum, and it was shown that in pure gauge theory these BC break the $`Z(3)`$ symmetry explicitely, which has important consequences for the high-temperature deconfinement phase transition. These conditions are also useful in numerical lattice simulations of this transition. When quarks are present, $`C`$-periodic BC break both chiral and flavour symmetries.

These boundary conditions are also especially important when topological properties are relevant in the system under consideration. This is the case of the lattice studies of confinement through monopole condensation. Recently the role of monopoles in connection with colour confinement has been evidentiated in SU(2) and SU(3) gluodynamics,? for which a disorder parameter based on the magnetic U(1) symmetry has been constructed and studied by Monte Carlo techniques. The disorder parameter is the vacuum expectation value (vev) of a disorder operator, which is an operator that creates a magnetic monopole in the gauge configuration. The definition of this disorder operator requires $`C`$-periodic BC in the time direction. For the pure gauge case, this means that the links at time $`t+N_t`$, where $`N_t`$ is the temporal extension of the lattice, are the complex conjugate of the links at time $`t`$. The effect on the simulation algorithm is a simple redefinition of the staples containing links that pierce the temporal boundary. The natural extension of the procedure used for the pure gauge case to full QCD requires the implementation of $`C`$-periodic BC in the presence of fermions. In particular, we will be concerned with the case of staggered fermions. $`C`$-periodic BC modify the fermionic matrix, and many proofs of properties used for the setup of standard simulation algorithms no longer hold.

$`C`$-periodic BC in the continuum are defined by the action of the charge conjugation operator $`C`$ on the fields.? However, lattice fermions are different from fermions in the continuum. In particular, $`C`$ is not a symmetry of the lattice action with staggered fermions. It also breaks translation invariance for finite lattice spacing. However, there is a discrete symmetry of the staggered fermion action,

$$S_f=\underset{i,\mu }{}\left[\frac{1}{2}\eta _{i,\mu }\left(\overline{\psi }_iU_{i,\mu }\psi _{i+\mu }\overline{\psi }_{i+\mu }U_{i,\mu }^{}\psi _i\right)+m\overline{\psi }_i\psi _i\right]$$
(1.1)
(here $`i`$ indicates the lattice point, $`U_{i,\mu }`$ is the SU(3) matrix associated with the link leaving the $`i`$-th lattice point in the $`\mu `$ direction, $`\psi _i`$ is the staggered fermion field at the point $`i`$, and $`\eta _{i,\mu }`$ the usual staggered fermion phase), which corresponds to charge conjugation in the gluon sector but which in the continuum limit also contains a flavour transformation.?,? We will call $`C^{}`$ this symmetry of the lattice action:

$${}_{}{}^{C^{}}U_{i,\mu }^{}=U_{i,\mu }^{},^C^{}\psi _i=ϵ_i\overline{\psi }_i^\mathrm{T},^C^{}\overline{\psi }_i=\psi _i^\mathrm{T}ϵ_i,$$
(1.2)
where $`ϵ_i=(1)^{x_i+y_i+z_i+t_i}`$, $`(x_i,y_i,z_i,t_i)`$ being the lattice coordinates of point $`i`$, the “T” represents the traspose operation, and we will use $`\psi ^{}\overline{\psi }^\mathrm{T}`$. Translation invariance implies that a BC must correspond to a symmetry of the action. $`C^{}`$-BC are defined as the boundary conditions corresponding to the symmetry (1.2):

$$\mathrm{\Phi }_{i+N}=^C^{}\mathrm{\Phi }_i,$$
(1.3)
where $`\mathrm{\Phi }`$ is a field, $`U_\mu `$ or $`\psi `$, and $`N`$ is the number of lattice points in the direction in which we use this boundary condition.

In the chiral limit the lattice action has a $`U(1)_\mathrm{E}U(1)_\mathrm{O}`$ chiral symmetry of independent rotations on $`(x+y+z+t)`$-even and -odd lattice points. $`C^{}`$-BC break explicitely this symmetry to $`U(1)_{\mathrm{E}=\mathrm{O}^{}}`$. Baryon number $`U(1)_{\mathrm{E}=\mathrm{O}}`$ is also broken explicitely down to $`Z(2)_{\mathrm{E}=\mathrm{O}}`$.?

In this paper we will have in mind the physical problem mentioned above: computation of the vev of a monopole creation operator in lattice QCD. This means that we will consider imposing $`C^{}`$-BC in the time direction, and periodic BC in the spatial directions. However this is just to fix the notation in what follows; the $`C^{}`$ conditions could in fact be assumed in any direction. The purpose of the paper is to show the theoretical framework to be used in a lattice simulation with these BC (section 2) and how usual algorithms? need to be modified (section 3). Our conclusions are summarized in section 4.

## 3 Mathematical description

Let us consider the partition function of lattice QCD with staggered fermions

$$Z=(𝒟U)(𝒟\psi 𝒟\overline{\psi })e^{S_gS_f},$$
(3.1)
where $`S_g(U)`$ is the Wilson action for the pure gauge sector, and $`S_f`$ is given by Eq. (1.1). The fermionic variables can be integrated out to give

$$Z=(𝒟U)e^{S_g(U)}detM(U),$$
(3.2)
where***We have absorbed the staggered phases by a redefinition of the link matrices $`U`$.

$$M(U)_{i,j}=m\delta _{i,j}+\underset{\mu }{}\frac{1}{2}(U_{i,\mu }\delta _{i,j\mu }U_{i\mu ,\mu }^{}\delta _{i,j+\mu })$$
(3.3)
is the fermionic matrix, and periodic BC are assumed. Using this matrix notation, we can write

$$S_f=\overline{\psi }M\psi =\frac{1}{2}\left[\left(\overline{\psi }M\psi \right)+\left(\overline{\psi }M\psi \right)^\mathrm{T}\right],$$
(3.4)
since $`S_f`$ is a number. We can use the new variable $`\mathrm{\Psi }`$ defined as the column vector formed by $`\psi `$ and $`\psi ^{}`$, so that the fermionic integral is

$$(𝒟\psi )(𝒟\overline{\psi })e^{S_f}=(𝒟\mathrm{\Psi })e^{\frac{1}{2}\mathrm{\Psi }^\mathrm{T}A\mathrm{\Psi }}=\mathrm{Pf}(A),$$
(3.5)
$$\mathrm{\Psi }^\mathrm{T}A\mathrm{\Psi }\left(\begin{array}{cc}\psi ^\mathrm{T}& \overline{\psi }\end{array}\right)A\left(\begin{array}{c}\psi \hfill \\ \psi ^{}\hfill \end{array}\right),$$
(3.6)
with

$$A=\left(\begin{array}{cc}0& M^\mathrm{T}\\ M& 0\end{array}\right),$$
(3.7)
and we have introduced the Pfaffian of the matrix $`A`$, $`\mathrm{Pf}(A)`$. It is well known that?

$$\mathrm{Pf}^2(A)=detA.$$
(3.8)
Now let us consider $`C^{}`$-BC in the time direction. Then, following Eqs. (1.2) and (1.3), $`\psi _{N_t}=ϵ_0\psi _0^{}`$, and $`\overline{\psi }_{N_t}=ϵ_0\psi _0^\mathrm{T}`$, where we have written in the subscript the temporal coordinate and omitted the spatial coordinates. Now Eq. (3.4) gives, for the terms connecting the slices $`N_t1`$ and $`N_t`$,

$`{\displaystyle \frac{1}{2}}\left(\overline{\psi }_{N_t1}U_{N_t1,t}\psi _0^{}+\psi _0^\mathrm{T}U_{N_t1,t}^{}\psi _{N_t1}\right)ϵ_0`$ (3.9)
$``$ $`{\displaystyle \frac{1}{2}}\left(\overline{\psi }_0U_{N_t1,t}^\mathrm{T}\psi _{N_t1}^{}+\psi _{N_t1}^\mathrm{T}U_{N_t1,t}^{}\psi _0\right)ϵ_0.`$
In this way the matrix $`A`$ of Eq. (3.7) is substituted by

$$A=\left(\begin{array}{cc}B& \stackrel{~}{M}^\mathrm{T}\\ \stackrel{~}{M}& B^{}\end{array}\right),$$
(3.10)
where $`B`$ satisfies the properties:

$$B^{}=B^{}B=B^\mathrm{T}$$
(3.11)
and $`\stackrel{~}{M}`$ is the fermionic matrix Eq. (3.3) apart from the terms connecting the slices $`N_t1`$ and $`N_t`$, which have gone to the matrices $`B`$ and $`B^{}`$.

Eq. (3.5) is still valid and we are interested in calculating $`\mathrm{Pf}(A)`$.

### 3.1 Pseudofermionic variables

The usual approach? to the simulation of theories with dynamical fermions is to rewrite the determinant of the fermionic matrix in Eq. (3.2) using that

$$det(M^{}M)𝒟\varphi ^{}𝒟\varphi \mathrm{exp}[\varphi ^{}(M^{}M)^1\varphi ],$$
(3.12)
where $`\varphi `$ is a complex bosonic field with the same quantum numbers as the Grassmann field. One introduces the matrix $`(M^{}M)^1`$, instead of $`M^1`$, so that the pseudofermionic fields can be generated using a simple heatbath method. Since $`detM`$ is a real number, $`det(M^{}M)=(detM)^2`$. So, actually, this corresponds to a double number of flavours with respect to the original theory. However, the matrix $`M^{}M`$ has two important properties: it has no matrix elements connecting even and odd lattice sites, and the determinants of its submatrices on the even and odd sites are equal.? Therefore one can avoid the redoubling of flavours by defining the pseudofermionic field only on even lattice sites.

A remarkable difference between the partition function with periodic BC and the partition function with $`C^{}`$-BC is that in the latter case the determinant of $`M`$, Eq. (3.2), has to be replaced by the Pfaffian of $`A`$, i.e. $`\pm \sqrt{detA}`$. Because of the square root, the usual trick of introducing pseudofermionic fields and rewriting this factor as the integral over these fields can only be applied if the number of continuum fermion flavours is such that the square root cancels. Moreover, in order to have a positive-definite integration measure, we need that the sign in front of this factor be $`+`$. Both conditions are satisfied if the number of continuum flavours is a multiple of eight. This generates a further unavoidable redoubling. Until Sect. 3.2, we will be concerned with the numerical simulation of $`detA`$, which is equivalent to simulating a double number of fermion flavours in the continuum limit (eight instead of four). In Sect. 3.2 we will discuss how to deal with the usual case of four staggered fermion flavours.

Since the case of a system with eight fermion flavours in the continuum limit and $`C^{}`$-BC and the case of a system with four fermion flavours in the continuum limit and periodic BC are similar, we would like to follow in the former case the standard procedure to obtain the determinant of the matrix (3.10). First, in the Appendix A it is shown that $`detA`$ is a real number. So we can also in this case use the matrix $`A^{}A`$ to introduce the pseudofermionic field, which will have now twice the number of components as in the usual case. Second, we will now see that the matrix $`A^{}A`$ does not connect even and odd lattice sites.

Using the form of the fermionic matrix (3.3) and the definitions of $`B`$ and $`\stackrel{~}{M}`$ given by Eqs. (3.9) and (3.10), we can split the blocks of the matrix $`A`$ in even and odd lattice sites:

$$A=\left(\begin{array}{cccc}0& \frac{1}{2}B_{\mathrm{eo}}& m& \frac{1}{2}D_{\mathrm{oe}}^\mathrm{T}\\ \frac{1}{2}B_{\mathrm{oe}}& 0& \frac{1}{2}D_{\mathrm{eo}}^\mathrm{T}& m\\ m& \frac{1}{2}D_{\mathrm{eo}}& 0& \frac{1}{2}B_{\mathrm{eo}}^{}\\ \frac{1}{2}D_{\mathrm{oe}}& m& \frac{1}{2}B_{\mathrm{oe}}^{}& 0\end{array}\right).$$
(3.13)
From Eqs. (3.3) and (3.11) it is easy to see that

$`D_{\mathrm{oe}}^{}=D_{\mathrm{eo}}`$ $`D_{\mathrm{oe}}^\mathrm{T}=D_{\mathrm{eo}}^{}`$ (3.14)
$`B_{\mathrm{oe}}^{}=B_{\mathrm{eo}}^{}`$ $`B_{\mathrm{oe}}^\mathrm{T}=B_{\mathrm{eo}}.`$ (3.15)
In the expression (3.13), the new blocks divide the matrix in rows and columns identified by the pairs $`[0,\mathrm{e}]`$, $`[0,\mathrm{o}]`$, $`[1,\mathrm{e}]`$ and $`[1,\mathrm{o}]`$, where the number indicates the blocks defined in Eq. (3.10) and the letter the even/odd subblock. Reorganizing the rows and columns to $`[0,\mathrm{e}]`$, $`[1,\mathrm{e}]`$, $`[0,\mathrm{o}]`$ and $`[1,\mathrm{o}]`$ (the determinant does not change), we rewrite the matrix A in the following form:

$$A=\left(\begin{array}{cc}\mu & \frac{1}{2}A_{\mathrm{eo}}\\ \frac{1}{2}A_{\mathrm{oe}}& \mu \end{array}\right),$$
(3.16)
where

$$\mu =\left(\begin{array}{cc}0& m\\ m& 0\end{array}\right),$$
(3.17)
$$A_{\mathrm{eo}}=\left(\begin{array}{cc}B_{\mathrm{eo}}& D_{\mathrm{eo}}^{}\\ D_{\mathrm{eo}}& B_{\mathrm{eo}}^{}\end{array}\right),$$
(3.18)
and $`A_{\mathrm{oe}}`$ has the same form as $`A_{\mathrm{eo}}`$ with the exchange $`\mathrm{e}\mathrm{o}`$. Eqs. (3.14) and (3.15) give

$$A_{\mathrm{oe}}^{}=A_{\mathrm{eo}}^{}.$$
(3.19)
From Eqs. (3.3) and (3.9) we obtain the explicit form of the matrices $`D_{\mathrm{eo}}`$ and $`B_{\mathrm{eo}}`$:

$$(D_{\mathrm{eo}})_{j,k}=\underset{\mu }{}(U_{j,\mu }\delta _{k,j+\mu }U_{j\mu ,\mu }^{}\delta _{k,j\mu }),$$
(3.20)
where $`j`$ is an even site and the link between $`j`$ and $`k`$ does not connect the $`t=N_t1`$ and $`t=0`$ time slices (otherwise $`(D_{\mathrm{eo}})_{j,k}=0`$), and

$$\begin{array}{cc}(B_{\mathrm{eo}})_{j,k}=U_{j,t}^{}\hfill & \hfill j=N_t1,k=0\\ (B_{\mathrm{eo}})_{j,k}=U_{k,t}^{}\hfill & \hfill j=0,k=N_t1\\ (B_{\mathrm{eo}})_{j,k}=0\hfill & \hfill \text{in every other case}\end{array}\}$$
(3.21)
where $`j=N_t1`$ or $`j=0`$ means that the temporal coordinate of site $`j`$ is $`N_t1`$ or zero, respectively. When $`j`$ is an odd site, $`(D_{\mathrm{oe}})_{j,k}`$ has the same expression as Eq. (3.20), and

$$\begin{array}{cc}(B_{\mathrm{oe}})_{j,k}=U_{j,t}^{}\hfill & \hfill j=N_t1,k=0\\ (B_{\mathrm{oe}})_{j,k}=U_{k,t}^{}\hfill & \hfill j=0,k=N_t1\\ (B_{\mathrm{oe}})_{j,k}=0\hfill & \hfill \text{in every other case}\end{array}\}$$
(3.22)
These expressions completely determine every element of the matrices $`A_{\mathrm{eo}}`$ and $`A_{\mathrm{oe}}`$.

We now compute

$$A^{}A=\left(\begin{array}{cc}\mu ^2+\frac{1}{4}A_{\mathrm{oe}}^{}A_{\mathrm{oe}}& \frac{1}{2}(\mu A_{\mathrm{eo}}A_{\mathrm{oe}}^{}\mu )\\ \frac{1}{2}(A_{\mathrm{eo}}^{}\mu \mu A_{\mathrm{oe}})& \mu ^2+\frac{1}{4}A_{\mathrm{eo}}^{}A_{\mathrm{eo}}\end{array}\right).$$
(3.23)
But using the properties (3.14) and (3.15) it is direct to see that the blocks out of the diagonal in Eq. (3.23) are zero. Then, using Eq. (3.19),

$$A^{}A=\left(\begin{array}{cc}\mu ^2\frac{1}{4}A_{\mathrm{eo}}^{}A_{\mathrm{oe}}& 0\\ 0& \mu ^2\frac{1}{4}A_{\mathrm{oe}}^{}A_{\mathrm{eo}}\end{array}\right)$$
(3.24)
and, as a result, we see that the matrix $`A^{}A`$ connects only lattice points of the same parity.

In order to use the same trick to avoid the flavour doubling produced by the introduction of $`A^{}A`$, that is, to define the pseudofermionic field on even sites only, we need to show that, also in this case, the determinant of the even and odd parts in the matrix (3.24) are equal. This is done in the following subsection.

### 3.2 Reducing the number of flavours: even-odd partitioning

Let us write $`KA^{}A`$. Then we have

$$K=\left(\begin{array}{cc}m^2I_2\frac{1}{4}A_{\mathrm{eo}}^{}A_{\mathrm{oe}}& 0\\ 0& m^2I_2\frac{1}{4}A_{\mathrm{oe}}^{}A_{\mathrm{eo}}\end{array}\right)\left(\begin{array}{cc}K_\mathrm{e}& 0\\ 0& K_\mathrm{o}\end{array}\right),$$
(3.25)
where $`I_2`$ is the $`2\times 2`$ identity matrix. We need to show that $`detK_\mathrm{e}=detK_\mathrm{o}`$. The strategy will be to show that both $`detK_\mathrm{e}`$ and $`detK_\mathrm{o}`$ are equal to $`detA`$. To this aim, let us consider the matrix $`A`$ written in the e–o partitioned form, Eq. (3.16),

$$A=\left(\begin{array}{cc}\mu & \frac{1}{2}A_{\mathrm{eo}}\\ \frac{1}{2}A_{\mathrm{oe}}& \mu \end{array}\right),$$
(3.26)
where

$$\mu =\left(\begin{array}{cc}0& m\\ m& 0\end{array}\right),$$
(3.27)
$$A_{\mathrm{eo}}=\left(\begin{array}{cc}B_{\mathrm{eo}}& D_{\mathrm{eo}}^{}\\ D_{\mathrm{eo}}& B_{\mathrm{eo}}^{}\end{array}\right),$$
(3.28)
and

$$A_{\mathrm{oe}}=\left(\begin{array}{cc}B_{\mathrm{oe}}& D_{\mathrm{oe}}^{}\\ D_{\mathrm{oe}}& B_{\mathrm{oe}}^{}\end{array}\right).$$
(3.29)
We will make use of a general property of the determinant of a square $`2N\times 2N`$ block matrix,

$$det\left(\begin{array}{cc}X& Y\\ W& Z\end{array}\right)=det(X)det(ZWX^1Y),$$
(3.30)
where $`X,Y,W,Z`$ are $`N\times N`$ square matrices and $`X`$ is invertible. Applying this equality to the computation of $`detA`$, we obtain

$$detA=det(\mu )det\left(\mu \frac{1}{4}A_{\mathrm{oe}}\mu ^1A_{\mathrm{eo}}\right)=det\left(\mu ^2\frac{1}{4}\mu A_{\mathrm{oe}}\mu ^1A_{\mathrm{eo}}\right).$$
(3.31)
Using the explicit form of $`\mu `$ and $`A_{\mathrm{oe}}`$ it is easily shown that

$$\mu A_{\mathrm{oe}}=A_{\mathrm{oe}}^{}\mu .$$
(3.32)
Inserting this equality in Eq. (3.31), it follows that

$$detA=det\left(\mu ^2+\frac{1}{4}A_{\mathrm{oe}}^{}A_{\mathrm{eo}}\right)=detK_\mathrm{o},$$
(3.33)
where the last equality is implied by the fact that $`N`$ is an even number.

We can compute again $`detA`$ using the same procedure after exchanging even and odd variables in $`A`$. Using a property analogous to that in Eq. (3.32), namely $`\mu A_{\mathrm{eo}}=A_{\mathrm{eo}}^{}\mu `$, it is then easily shown that

$$detA=det\left(\mu ^2+\frac{1}{4}A_{\mathrm{eo}}^{}A_{\mathrm{oe}}\right)=detK_\mathrm{e}.$$
(3.34)
This is a proof that $`detK_\mathrm{e}=detK_\mathrm{o}`$.

We have shown that $`detK_\mathrm{e}detK_\mathrm{o}=(detA)^2`$. On the other hand it is also true that $`detK_\mathrm{e}detK_\mathrm{o}=det(A^{}A)=detA^{}detA`$, therefore we have also obtained a proof that $`detA=detA^{}`$ alternative to that given in Appendix A. Moreover, we have that

$$detA=detK_\mathrm{e}=det\left(m^2I_2+\frac{1}{4}A_{\mathrm{oe}}^{}A_{\mathrm{oe}}\right)>0,$$
(3.35)
since both $`m^2I_2`$ and $`A_{\mathrm{oe}}^{}A_{\mathrm{oe}}`$ are positive definite matrices. Therefore $`detA`$ is a positive number.

We notice that this method can be easily applied also to the standard case, thus providing a proof alternative to those presented in Ref. 8. This is shown in detail in Appendix B.

## 4 Hybrid Monte Carlo Implementation

We will show now how the standard Hybrid Monte Carlo (HMC) algorithm? needs to be modified to incorporate $`C^{}`$-BC. In this algorithm one introduces fictitious momenta, conjugate variables of the links, as dynamical variables, and makes fields evolve with a mixed dynamics, in which deterministic and stochastic steps are alternated in a prescribed way. In the deterministic part of the algorithm, the system follows the equations of motion derived from the Hamiltonian of the (4+1)-dimensional system. These equations give the evolution of the fields in the fictitious time $`\tau `$. The equation of motion for the matrix $`U_{j,\mu }`$ is

$$\dot{U}_{j,\mu }=iH_{j,\mu }U_{j,\mu },$$
(4.1)
where the conjugate momentum $`H_{j,\mu }`$ is a traceless Hermitian matrix, and $`\dot{U}`$ is the derivative of $`U`$ with respect to $`\tau `$. The equations of motion for the momenta $`H`$ are obtained by imposing that the Hamiltonian be constant. The integration of these equations of motion is carried out numerically after discretization of $`\tau `$. Usually the temporal step is of order $`10^2`$$`10^3`$ and the configuration space is sampled with a total length of the trajectory of order 1. The stochastic part of the algorithm consists in the generation of new momenta and new pseudofermionic variables according to their probability distributions at the beginning of each trajectory. Moreover, at the end of each trajectory a Metropolis accept-reject step is performed, which makes the algorithm exact.

### 4.1 HMC algorithm with $`C^{}`$ boundary conditions

Once introduced the pseudofermionic fields, defined only on even sites, and the auxiliary momenta fields, the partition function of the system is

$$Z=(𝒟U𝒟\varphi ^{}𝒟\varphi 𝒟H)e^{},$$
(4.2)
with

$$=\frac{1}{2}\underset{j,\mu }{}\mathrm{tr}H_{j,\mu }^2+S_g+\varphi ^{}(A^{}A)^1\varphi .$$
(4.3)
To obtain an equation of motion for $`H`$ we require that $``$ be a constant of motion, that is, $`\dot{}=0`$. The tricky part in this differentiation is in the fermionic term. The derivative of $`_f=\varphi ^{}(A^{}A)^1\varphi `$ is

$`\dot{}_f`$ $`=`$ $`{\displaystyle \underset{j,\mu }{}}\varphi ^{}(A^{}A)^1i[{\displaystyle \frac{A^{}}{U_{j,\mu }}}H_{j,\mu }U_{j,\mu }A+{\displaystyle \frac{A^{}}{U_{j,\mu }^\mathrm{T}}}U_{j,\mu }^\mathrm{T}H_{j,\mu }^\mathrm{T}A+`$ (4.4)
$`A^{}{\displaystyle \frac{A}{U_{j,\mu }}}H_{j,\mu }U_{j,\mu }+A^{}{\displaystyle \frac{A}{U_{j,\mu }^\mathrm{T}}}U_{j,\mu }^\mathrm{T}H_{j,\mu }^\mathrm{T}H_{j,\mu }^\mathrm{T}U_{j,\mu }^{}{\displaystyle \frac{A^{}}{U_{j,\mu }^{}}}A`$
$`U_{j,\mu }^{}H_{j,\mu }{\displaystyle \frac{A^{}}{U_{j,\mu }^{}}}AA^{}H_{j,\mu }^\mathrm{T}U_{j,\mu }^{}{\displaystyle \frac{A}{U_{j,\mu }^{}}}A^{}U_{j,\mu }^{}H_{j,\mu }{\displaystyle \frac{A}{U_{j,\mu }^{}}}]`$
$`\times (A^{}A)^1\varphi ,`$
where we have made use of Eq. (4.1) and the fact that $`H`$ is Hermitian. We have also used that $`A`$ is linear in $`U`$, $`U^\mathrm{T}`$, $`U^{}`$ and $`U^{}`$, and then the partial derivatives in the previous expression commute with $`H`$ and $`U`$.

It is convenient to introduce the operator

$$P_{ij}=X_iX_j^{},$$
(4.5)
where

$$X=(A^{}A)^1\varphi .$$
(4.6)
It is easily seen that $`P`$ is Hermitian:

$$(P^{})_{ij}=(P_{ji})^{}=X_iX_j^{}=P_{ij},$$
(4.7)
and, being $`\varphi `$ defined only on even sites, $`P_{ij}`$ is taken to be zero unless $`i`$ and $`j`$ are both even sites. On the other hand, we have that

$`\left({\displaystyle \frac{A^{}}{U}}\right)^{}={\displaystyle \frac{A}{U^{}}}`$ $`\left({\displaystyle \frac{A}{U}}\right)^{}={\displaystyle \frac{A^{}}{U^{}}}`$
$`\left({\displaystyle \frac{A^{}}{U^\mathrm{T}}}\right)^{}={\displaystyle \frac{A}{U^{}}}`$ $`\left({\displaystyle \frac{A}{U^\mathrm{T}}}\right)^{}={\displaystyle \frac{A^{}}{U^{}}}`$ (4.8)
and then we can write

$`\dot{}_f`$ $`=`$ $`{\displaystyle \underset{j,\mu }{}}\mathrm{tr}[iH_{j,\mu }(U_{j,\mu }AP{\displaystyle \frac{A^{}}{U_{j,\mu }}}+U_{j,\mu }PA^{}{\displaystyle \frac{A}{U_{j,\mu }}})+`$ (4.9)
$`iH_{j,\mu }^\mathrm{T}(AP{\displaystyle \frac{A^{}}{U_{j,\mu }^\mathrm{T}}}U_{j,\mu }^\mathrm{T}+PA^{}{\displaystyle \frac{A}{U_{j,\mu }^\mathrm{T}}}U_{j,\mu }^\mathrm{T})+\mathrm{H}.\mathrm{c}.],`$
where H.c. means the Hermitian conjugate and we have used the cyclic property of the trace operation.

The trace in Eq. (4.9) is taken over both color and site indices. Taking this into account and using the cyclic property again, we have

$`\dot{}_f`$ $`=`$ $`{\displaystyle \underset{j,\mu }{}}\mathrm{tr}(H_{j,\mu }\{i[U_{j,\mu }\left(AP{\displaystyle \frac{A^{}}{U_{j,\mu }}}\right)+U_{j,\mu }\left(PA^{}{\displaystyle \frac{A}{U_{j,\mu }}}\right)+`$ (4.10)
$`U_{j,\mu }\left(AP{\displaystyle \frac{A^{}}{U_{j,\mu }^\mathrm{T}}}\right)^\mathrm{T}+U_{j,\mu }\left(PA^{}{\displaystyle \frac{A}{U_{j,\mu }^\mathrm{T}}}\right)^\mathrm{T}]+\mathrm{H}.\mathrm{c}.\}).`$
The following step is to calculate the derivatives appearing in Eq. (4.10). From Eqs. (3.16) and (3.19), we have that

$$A^{}=A^{}.$$
(4.11)
We write again $`A`$ and $`A^{}`$ in terms of subblocks, identifying the rows and columns with the pairs $`[0,\mathrm{e}]`$, $`[1,\mathrm{e}]`$, $`[0,\mathrm{o}]`$ and $`[1,\mathrm{o}]`$:

$$A=\left(\begin{array}{cccc}0& m& \frac{1}{2}B_{\mathrm{eo}}& \frac{1}{2}D_{\mathrm{eo}}^{}\\ m& 0& \frac{1}{2}D_{\mathrm{eo}}& \frac{1}{2}B_{\mathrm{eo}}^{}\\ \frac{1}{2}B_{\mathrm{oe}}& \frac{1}{2}D_{\mathrm{oe}}^{}& 0& m\\ \frac{1}{2}D_{\mathrm{oe}}& \frac{1}{2}B_{\mathrm{oe}}^{}& m& 0\end{array}\right),$$
(4.12)
$$A^{}=\left(\begin{array}{cccc}0& m& \frac{1}{2}B_{\mathrm{eo}}^{}& \frac{1}{2}D_{\mathrm{eo}}\\ m& 0& \frac{1}{2}D_{\mathrm{eo}}^{}& \frac{1}{2}B_{\mathrm{eo}}\\ \frac{1}{2}B_{\mathrm{oe}}^{}& \frac{1}{2}D_{\mathrm{oe}}& 0& m\\ \frac{1}{2}D_{\mathrm{oe}}^{}& \frac{1}{2}B_{\mathrm{oe}}& m& 0\end{array}\right).$$
(4.13)
Once written everything in terms of $`B`$ and $`D`$ we can easily see where $`U`$ and $`U^\mathrm{T}`$ appear, from Eqs. (3.20)–(3.22), and compute the derivatives. The result is given in Table 4.1.

Now it is easy to calculate the different terms in Eq. (4.10). Let us do, as an example, the first one. The trace over lattice site indices affects only the expressions between parentheses. From Table 4.1, we see that, for $`j`$ even, and $`jN_t1`$ or $`\mu t`$, the first one gives

$`\mathrm{tr}\left(AP{\displaystyle \frac{A^{}}{U_{j,\mu }}}\right)`$ $`=`$ $`A_{k_2,k}X_kX_{k_1}^{}\left({\displaystyle \frac{A^{}}{U_{j,\mu }}}\right)_{k_1,k_2}={\displaystyle \frac{1}{2}}A_{[1,\mathrm{o}]j+\mu ,k}X_kX_{[0,\mathrm{e}]j}^{}=`$ (4.14)
$`\left({\displaystyle \frac{1}{4}}(D_{\mathrm{oe}})_{j+\mu ,k}X_k^0+{\displaystyle \frac{1}{4}}(B_{\mathrm{oe}}^{})_{j+\mu ,k}X_k^1\right)(X^{})_j^0=`$
$`{\displaystyle \frac{1}{4}}((A_{\mathrm{oe}})_{j+\mu ,k}X_k)^1(X^{})_j^0={\displaystyle \frac{1}{4}}(A_{\mathrm{oe}}X)_{j+\mu }^1(X^{})_j^0,`$
where we have used Eqs. (4.12) and (3.18), and called $`X_k^0`$ the first three components of the vector $`X_k`$, and $`X_k^1`$ the second three components. Notice that there is not a mass term because it is diagonal, which would imply that $`k=j+\mu `$, that is, $`k`$ would be odd, and then $`X_k=0`$. For $`j`$ even and $`j=N_t1`$, $`\mu =t`$, an analogous calculation gives

$$\mathrm{tr}\left(AP\frac{A^{}}{U_{j,\mu }}\right)=\frac{1}{4}(A_{\mathrm{oe}}X)_{j+\mu }^0(X^{})_j^0.$$
(4.15)
This term does not contribute when $`j`$ is odd, because in this case $`(X^{})_j^0`$ is zero ($`X`$ is defined on even sites only).

Computing every other contribution in Eq. (4.10), the final result can be written

$$\dot{}_f=\underset{j,\mu }{}\mathrm{tr}\left\{H_{j,\mu }\left(iF_{j,\mu }iF_{j,\mu }^{}\right)\right\},$$
(4.16)
where $`F_{j,\mu }`$ takes a different form depending on $`j`$ and $`\mu `$.

For $`j`$ even and $`(j,\mu )(N_t1,t)`$,

$$F_{j,\mu }=U_{j,\mu }\left(\frac{1}{4}(A_{\mathrm{oe}}X)_{j+\mu }^1(X^{})_j^0\left(\frac{1}{4}X_j^1((A_{\mathrm{oe}}X)^{})_{j+\mu }^0\right)^\mathrm{T}\right);$$
(4.17)
for $`j`$ even and $`(j,\mu )=(N_t1,t)`$,

$$F_{j,\mu }=U_{j,\mu }\left(\frac{1}{4}(A_{\mathrm{oe}}X)_{j+\mu }^0(X^{})_j^0+\left(\frac{1}{4}X_j^1((A_{\mathrm{oe}}X)^{})_{j+\mu }^1\right)^\mathrm{T}\right);$$
(4.18)
for $`j`$ odd and $`(j,\mu )(N_t1,t)`$,

$$F_{j,\mu }=U_{j,\mu }\left(\frac{1}{4}X_{j+\mu }^0((A_{\mathrm{oe}}X)^{})_j^1+\left(\frac{1}{4}(A_{\mathrm{oe}}X)_j^0(X^{})_{j+\mu }^1\right)^\mathrm{T}\right);$$
(4.19)
for $`j`$ odd and $`(j,\mu )=(N_t1,t)`$,

$$F_{j,\mu }=U_{j,\mu }\left(\frac{1}{4}X_{j+\mu }^1((A_{\mathrm{oe}}X)^{})_j^1\left(\frac{1}{4}(A_{\mathrm{oe}}X)_j^0(X^{})_{j+\mu }^0\right)^\mathrm{T}\right).$$
(4.20)
Since we determined completely the elements of the matrices $`A_{\mathrm{eo}}`$ and $`A_{\mathrm{oe}}`$ in section 2.1, these expressions allow us to obtain $`F_{j,\mu }`$ as a function of the links $`U`$ for every $`j`$ and $`\mu `$.

Coming back to Eq. (4.3), the condition $`\dot{}=0`$ gives

$$\dot{}=0=\underset{j,\mu }{}\mathrm{tr}(\dot{H}_{j,\mu }H_{j,\mu }+\frac{\beta }{6}(iH_{j,\mu }U_{j,\mu }V_{j,\mu }+\mathrm{H}.\mathrm{c}.)(iH_{j,\mu }F_{j,\mu }+\mathrm{H}.\mathrm{c}.)),$$
(4.21)
where $`V_{j,\mu }`$ is the sum of staples, or products of the other three matrices in the plaquettes containing $`U_{j,\mu }`$, and arises from the differentiation of the Wilson action $`S_g`$. Of course, the staples at the border of the lattice contain links defined by the boundary conditions ($`C^{}`$ or periodic, depending on the direction). The final solution for $`\dot{H}_{j,\mu }`$ is

$$i\dot{H}_{j,\mu }=\left[\frac{\beta }{3}U_{j,\mu }V_{j,\mu }2F_{j,\mu }\right]_{\mathrm{T}A},$$
(4.22)
where the subscript TA indicates the traceless anti-Hermitian part of the matrix:

$$Q_{\mathrm{T}A}=\frac{1}{2}(QQ^{})\frac{1}{6}\mathrm{tr}(QQ^{}).$$
(4.23)
### 4.2 Reducing the number of flavours: the Hybrid algorithm

Because of Eq. (3.8), the HMC algorithm that we have just described simulates eight fermion flavours in the continuum. In order to come back to four flavours, we note that

$$\mathrm{Pf}(A)=\pm (detA)^{1/2}$$
(4.24)
(we saw in Eq. (3.35) that $`detA`$ is a positive number).

One can simulate $`(detA)^{1/2}`$ by reverting to an approximate algorithm. A suitable choice is the $`R`$ algorithm,? in which discretization errors in the molecular dynamics part are of $`𝒪(\mathrm{\Delta }\tau ^2)`$.

On a lattice closed with $`C^{}`$-BC, the full QCD action in the presence of $`N_f`$ families of degenerate continuum fermions can be written as follows:

$$Z=(𝒟U)[det(A^{}A)_\mathrm{e}]^{N_f/8}e^{S_g(U)}=(𝒟U)e^{S_g(U)+\frac{N_f}{8}\mathrm{tr}\mathrm{log}(A^{}A)_\mathrm{e}},$$
(4.25)
where $`(A^{}A)_\mathrm{e}`$ is the restriction of $`A^{}A`$ to the even lattice sites.

Since the implementation of the $`R`$ algorithm in the present case reduces to simple adaptation of a standard technique to the system described by the equation of motion obtained in the previous subsection, we do not elaborate further on this point. However, we still have the “sign problem” of Eq. (4.24). The sign of the Pfaffian can be included by reweighting the expectation values according to

$$𝒪=\frac{𝒪\mathrm{signPf}(A)_+}{\mathrm{signPf}(A)_+},$$
(4.26)
where $`\mathrm{}_+`$ means that the expectation values have been obtained simulating $`+(detA)^{1/2}`$. One then needs to monitorize the sign of the Pfaffian depending on the gauge configurations. Some techniques to do that are explained in Ref. 9.

## 5 Conclusions

$`C^{}`$ boundary conditions are interesting to study some spontaneous symmetry breaking aspects of QCD. They are relevant when one analyses confinement through monopole condensation. We have shown in this work how these boundary conditions can be implemented to carry out a lattice simulation of full QCD with staggered fermions. We have proved that the common even-odd trick used to avoid the fermion redoubling produced by the introduction of the pseudofermionic field can be applied also to this case. However, there is an additional redoubling which forces to work with a minimum number of eight flavours with the usual Hybrid Monte Carlo algorithm, which we have adapted to this case. An alternative to avoid that is to consider a non-exact algorithm, which can be applied to any number of flavours.

These algorithms have been implemented and are presently running on an APE Quadrics machine to explore monopole condensation in full QCD.

Acknowledgments

This work has been partially supported by EU TMR program ERBFMRX-CT97-0122, Italian MURST and PPARC Grant PPA/G/0/1998/00567. We thank valuable comments from Isabel Campos, Simon Hands, Pilar Hernández and Giampiero Paffuti.

References

Appendix A

We will show that $`detA`$ is a real number, or, equivalently, that $`detA^{}=detA`$. This property is required in order that one can use the matrix $`A^{}A`$ to introduce the pseudofermionic field.

Let us consider $`A`$ in the form

$$A=\left(\begin{array}{cc}B& mD^\mathrm{T}\\ m+D& B^{}\end{array}\right),$$
(A.1)
where $`D_{i,j}=\stackrel{~}{M}_{i,j}m\delta _{i,j}`$.

Using the general property of the determinant of a square block matrix reported in Eq. (3.30) and exchangingThis corresponds to an even number of row exchanges, since $`A`$ has dimension $`2N`$ and $`N`$, being the dimension of the fermion matrix with periodic boundary conditions, is always an even number. Therefore $`detA`$ does not change under this operation. the two rows of blocks of $`A`$ we can rewrite

$`detA`$ $`=`$ $`det\left(\begin{array}{cc}m+D& B^{}\\ B& mD^\mathrm{T}\end{array}\right)`$ (A.4)
$`=`$ $`det(m+D)det(mD^\mathrm{T}B(D+m)^1B^{}).`$ (A.5)
Remembering that $`D^{}=D`$, we can write $`A^{}`$ in the form

$$A^{}=\left(\begin{array}{cc}B^{}& mD\\ m+D^\mathrm{T}& B\end{array}\right),$$
(A.6)
so that, exchanging the two columns of blocks in $`A^{}`$ and using again Eq. (3.30), we can write

$`detA^{}`$ $`=`$ $`det\left(\begin{array}{cc}mD& B^{}\\ B& m+D^\mathrm{T}\end{array}\right)`$ (A.9)
$`=`$ $`det(mD)det(m+D^\mathrm{T}B(D+m)^1B^{}).`$ (A.10)
After extracting a factor $`(1)^Nm^{2N}=m^{2N}`$ from both $`detA`$ and $`detA^{}`$ and defining $`\alpha =1/m`$, we can write

$`detA`$ $`=`$ $`m^{2N}det(1+\alpha D)det(1+\alpha D^\mathrm{T}+\alpha ^2B(1+\alpha D)^1B^{})`$
$`detA^{}`$ $`=`$ $`m^{2N}det(1\alpha D)det(1\alpha D^\mathrm{T}+\alpha ^2B(1\alpha D)^1B^{})`$ (A.11)
It clearly appears from Eq. (A.11) that $`detA^{}`$ is obtained from $`detA`$ by changing the sign of $`\alpha `$. Therefore, in order to show that $`detA=detA^{}`$, it is sufficient to show that $`detA`$ is an even function of $`\alpha `$.

Let us consider $`det(1+\alpha D)`$ at first. We can expand the determinant as follows:

$`det(1+\alpha D)`$ $`=`$ $`\mathrm{exp}\left(\mathrm{tr}\mathrm{ln}(1+\alpha D)\right)`$ (A.12)
$`=`$ $`\mathrm{exp}\left(\mathrm{tr}{\displaystyle \underset{k=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{(1)^k}{k}}(\alpha D)^k\right).`$
It is easy to see that the trace of the product of an odd number of $`D`$ matrices is zero. Indeed $`D`$ only connects nearest neighbour lattice sites, so it is not possible to connect a site to itself using the product of an odd number of $`D`$ matrices. Therefore only even powers of $`\alpha `$ appear in the expansion in Eq. (A.12) and $`det(1+\alpha D)`$ is an even function of $`\alpha `$.

Let us consider now $`det(1+\alpha D^\mathrm{T}+\alpha ^2B(1+\alpha D)^1B^{})`$, which we rewrite as $`det(1+P(\alpha ))`$, where

$$P(\alpha )=\alpha D^\mathrm{T}+\alpha ^2B\left(\underset{k=0}{\overset{\mathrm{}}{}}(1)^k\alpha ^kD^k\right)B^{}.$$
(A.13)
We notice that the matrix $`P(\alpha )`$ is expressed as series expansion in $`\alpha `$, where the coefficient of the $`k`$-th term is a homogeneous polynomial of degree $`k`$ in the matrices $`B,B^{},D`$ and $`D^\mathrm{T}`$. Therefore, expanding again the determinant as

$`det(1+P(\alpha ))`$ $`=`$ $`\mathrm{exp}\left(\mathrm{tr}\mathrm{ln}(1+P(\alpha ))\right)`$ (A.14)
$`=`$ $`\mathrm{exp}\left(\mathrm{tr}{\displaystyle \underset{k=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{(1)^k}{k}}P(\alpha )^k\right),`$
we see that $`det(1+P(\alpha ))`$ can be expanded as a power series in $`\alpha `$ and that the coefficient of the $`k`$-th term is the trace of a homogeneous polynomial of degree $`k`$ in the matrices $`B,B^{},D`$ and $`D^\mathrm{T}`$: since these matrices only connect nearest neighbour sites, the trace is zero for $`k`$ odd. Therefore also in this case the determinant is an even function of $`\alpha `$.

We conclude that $`detA`$, being the product of even functions of $`\alpha `$, is also an even function of $`\alpha `$, and therefore, from Eq. (A.11), $`detA=detA^{}`$.

Appendix B

We will give here a proof, alternative to that presented in Ref. 8, of the equality of the determinants of the submatrices on even and odd sites of $`M^{}M`$ in the case with standard boundary conditions.

In the standard case the fermion matrix $`M`$, defined in Eq. (3.3), has the following form

$$M=\left(\begin{array}{cc}m& \frac{1}{2}D_{\mathrm{eo}}\\ \frac{1}{2}D_{\mathrm{oe}}& m\end{array}\right),$$
(B.1)
and it is easily shown that $`D_{\mathrm{eo}}^{}=D_{\mathrm{oe}}`$. Using this last property it follows that

$$M^{}M=\left(\begin{array}{cc}m^2\frac{1}{4}D_{\mathrm{eo}}D_{\mathrm{oe}}& 0\\ 0& m^2\frac{1}{4}D_{\mathrm{oe}}D_{\mathrm{eo}}\end{array}\right).$$
(B.2)
Using the decomposition of the determinant given in Eq. (3.30), it is easy to show that

$$detM=det\left(m^2\frac{1}{4}D_{\mathrm{eo}}D_{\mathrm{oe}}\right).$$
(B.3)
On the other hand, if we exchange even and odd variables in $`M`$ before applying Eq. (3.30), we obtain

$$detM=det\left(m^2\frac{1}{4}D_{\mathrm{oe}}D_{\mathrm{eo}}\right),$$
(B.4)
and therefore

$$det\left(m^2\frac{1}{4}D_{\mathrm{oe}}D_{\mathrm{eo}}\right)=det\left(m^2\frac{1}{4}D_{\mathrm{eo}}D_{\mathrm{oe}}\right).$$
(B.5)


text:
# Swelling-Collapse Transition of Self-Attracting Walks

## Abstract

We study the structural properties of self-attracting walks in $`d`$ dimensions using scaling arguments and Monte Carlo simulations. We find evidence for a transition analogous to the $`\mathrm{\Theta }`$ transition of polymers. Above a critical attractive interaction $`u_\mathrm{c}`$, the walk collapses and the exponents $`\nu `$ and $`k`$, characterising the scaling with time $`t`$ of the mean square end-to-end distance $`R^2t^{2\nu }`$ and the average number of visited sites $`St^k`$, are universal and given by $`\nu =1/(d+1)`$ and $`k=d/(d+1)`$. Below $`u_\mathrm{c}`$, the walk swells and the exponents are as with no interaction, i.e. $`\nu =1/2`$ for all $`d`$, $`k=1/2`$ for $`d=1`$ and $`k=1`$ for $`d2`$. At $`u_\mathrm{c}`$, the exponents are found to be in a different universality class.

In recent years different models of random walks with memory or interaction have been studied. They can be divided in static  and dynamic  models, for an overview we refer to the papers of Duxbury and Queiroz  and Oettinger . Most efforts concentrated on models with repulsive interactions, in particular self-avoiding walks (SAW), which have been found useful for investigating polymers in dilute solution. When an attraction term $`\mathrm{exp}(A/T)`$, $`A<0`$, is included, the SAW model reveals a swelling-collapse transition at the ‘$`\mathrm{\Theta }`$ point’ $`T=\mathrm{\Theta }`$ . In contrast, the likewise challenging case of random walks with a similar attractive interaction, but without repulsion, has been less understood. This problem was solved only for one dimension, while in higher dimensions the results are highly controversial. Our numerical and analytical study of attractive random walks suggests that there exists a swelling-collapse transition, too, that is analogous to the $`\mathrm{\Theta }`$ transition in polymers.

We focus on the dynamic model of self-attracting walks (SATW, where a random walker jumps with probability $`p\mathrm{exp}(nu)`$  to a nearest neighbor site, with $`n=1`$ for already visited sites and $`n=0`$ for not visited sites. The interaction parameter $`u`$ is equivalent to $`A/T`$ for linear polymers. For $`u>0`$, the walk is attracted to its own trajectory . The structural behavior of the walk can be characterized by the mean square end-to-end distance $`R^2(t)`$ and the average number of visited sites $`S(t)`$. It is expected that these quantities scale with time $`t`$ as

$`R^2(t)t^{2\nu }(1\text{a})`$ (1)
and (2)
$`S(t)t^k.(3\text{b})`$ (3)
Earlier analyses for the SATW in two and three dimensions were not conclusive and the numerical data have been controversially interpreted . While Sapozhnikov  considered the possibility of the existence of a critical attraction $`u_\mathrm{c}`$ (but his numerical results were not conclusive), Lee  and Reis  argue strongly against the existence of $`u_\mathrm{c}`$, since they find $`\nu `$ and $`k`$ continuously decreasing with $`u`$.

In this Letter we present scaling arguments and extensive numerical simulations for $`R^2`$ and $`S`$ that strongly suggest the existence of a critical attraction $`u_\mathrm{c}`$ in $`d2`$, with three different universality classes for $`u>u_\mathrm{c}`$, $`u<u_\mathrm{c}`$ and $`u=u_\mathrm{c}`$. Below $`u_\mathrm{c}`$, the SATW is in the universality class of random walks, with $`\nu =1/2`$ and $`k=1`$. Above $`u_\mathrm{c}`$, the SATW collapses and the exponents change to $`\nu =1/(d+1)`$ and $`k=d/(d+1)`$. At the critical point, the exponents are $`\nu _\mathrm{c}=0.40\pm 0.01`$ and $`k_\mathrm{c}=0.80\pm 0.01`$ in $`d=2`$ and $`\nu _\mathrm{c}=0.32\pm 0.01`$ and $`k_\mathrm{c}=0.91\pm 0.03`$ in $`d=3`$ . The existence of $`u_\mathrm{c}`$ is in striking similarity to the ‘$`\mathrm{\Theta }`$ point’ phenomenon of linear polymers  where three different universality classes for $`T>\mathrm{\Theta }`$, $`T=\mathrm{\Theta }`$ and $`T<\mathrm{\Theta }`$ exist.

We used Monte-Carlo simulations to study $`R^2(t)`$ and $`S(t)`$. Figure 1 shows representative results of $`R^2(t)`$, for several values of $`u`$ in $`d=3`$. For large values of $`u`$, the curves bend down towards the slope of $`2\nu 0.5`$, while for small values of $`u`$, the curves bend up towards the slope of $`2\nu 1`$. At some intermediate critical value $`u_\mathrm{c}1.9`$, the slope is approximately $`2\nu _\mathrm{c}0.64`$. The mean number of visited sites $`S(t)`$ shows similar behavior, with $`k1`$ below $`u_\mathrm{c}`$, $`k_\mathrm{c}0.91`$ at $`u_\mathrm{c}`$, and $`k0.75`$ above $`u_\mathrm{c}`$. Figure 2a summarizes the asymptotic exponents $`\nu `$ and $`k`$ as a function of $`u`$ in $`d=3`$. We obtained similar results in $`d=2`$, the asymptotic values of $`\nu `$ and $`k`$ are presented in Fig. 2b. In Table I the values of the exponents are summarized and compared with the analogous known exponents for the $`\mathrm{\Theta }`$ transition in linear polymers.

In the following we present analytical arguments for the exponents above criticaltity, which can explain our numerical findings. We assume that for sufficiently strong attraction $`u>u_\mathrm{c}`$ the grown clusters are compact, so that the average number of visited sites scales with the rms displacement $`R(t)R^2(t)^{1/2}`$ as

$$S(t)R(t)^d,uu_\mathrm{c}.$$
(5)
Comparing Eq. (1) and Eq. (5) yields

$$k=\nu d,uu_\mathrm{c}.$$
(6)
For sufficiently strong attraction it takes a very long time for the walker to jump to an unvisited site. Before doing this, the walker diffuses around on the visited sites, being located with equal probability on any of the cluster sites. Hence the mean cluster growth rate is proportional to the ratio between the number of boundary sites and the total number of the cluster sites :

$$\frac{\mathrm{d}S}{\mathrm{d}t}\frac{R^{d1}}{R^d}t^\nu .$$
(7)
Thus $`St^{\nu +1}`$. Combining this result with Eq. (1b) and (6), we obtain

$`\nu ={\displaystyle \frac{1}{d+1}}(8\text{a})`$ (8)
and (9)
$`k={\displaystyle \frac{d}{d+1}}(10\text{b})`$ (10)
for $`uu_\mathrm{c}`$.

Because of universality we assume that these results, which are in agreement with the exact values $`\nu =1/2`$ and $`k=1/2`$ in $`d=1`$  and are supported by our extensive Monte Carlo simulations in $`d=2`$ and $`d=3`$, are valid for all $`u>u_\mathrm{c}`$. Indeed, Fig. 2 suggests that the predictions for $`u>u_\mathrm{c}`$ (Eq. (8)) are approached asymptotically. We like to note that in $`d=2`$ the relation $`k=\nu d`$ also holds for $`uu_\mathrm{c}`$, while in $`d=3`$ the numerical results yield $`k<\nu d`$ for $`u<u_\mathrm{c}`$. Since the mass of the generated clusters scales like $`MSR^{k/\nu }`$, $`k/\nu `$ corresponds to the fractal dimension $`d_\mathrm{f}`$ of the cluster. In $`d=2`$ the clusters are compact for all $`u`$ as $`k/\nu =d_\mathrm{f}=d`$. In $`d=3`$ they are compact for $`u>u_\mathrm{c}`$, while for $`u<u_\mathrm{c}`$, the fractal dimension of clusters generated by simple random walks $`d_\mathrm{f}=2<d`$ is obtained. At the criticality, we find $`d_\mathrm{f}=2.84\pm 0.25`$, but we cannot rule out the possibility that $`d_\mathrm{f}=d`$.

To understand the behavior in the critical regime we suggest the following scaling approach. Guided by Fig. 1, we assume that there exists a crossover time $`t_\xi `$ below which the exponent $`\nu `$ is close to $`\nu _\mathrm{c}`$ and above which $`\nu `$ approaches $`1/2`$ for $`u<u_\mathrm{c}`$ and $`1/(d+1)`$ for $`u>u_\mathrm{c}`$. This suggests the following scaling relations:

$`R(t)t^{\nu _\mathrm{c}}f_\pm (t/t_\xi )(12\text{a})`$ (12)
and (13)
$`S(t)t^{k_\mathrm{c}}g_\pm (t/t_\xi ),(14\text{b})`$ (14)
where (15)
$`t_\xi =|uu_\mathrm{c}|^\alpha .(16\text{c})`$ (16)
The plus sign refers to $`u>u_\mathrm{c}`$, the minus sign to $`u<u_\mathrm{c}`$, and the exponent $`\alpha `$ has to be determined numerically. As $`t_\xi `$ is assumed to be the only relevant time scale, the scaling functions bridge the short time and the long time regime. To match both regimes, we require that $`f_\pm (x)=\mathrm{const}`$ for $`x1`$ ($`tt_\xi `$), and $`f_+(x)x^{1/(d+1)\nu _\mathrm{c}}`$, $`f_{}(x)x^{1/2\nu _\mathrm{c}}`$ for $`x1`$. Analogous results are expected for $`g_\pm (x)`$, with $`g_\pm (x)=\mathrm{const}`$ for $`x1`$, and $`g_+(x)x^{d/(d+1)k_\mathrm{c}}`$, $`g_{}(x)x^{1k_\mathrm{c}}`$ for $`x1`$.

To test the scaling theory and to determine the exponent $`\alpha `$ we plotted $`R^2(t)/t_\xi ^{2\nu _\mathrm{c}}`$ and $`S(t)/t_\xi ^{k_\mathrm{c}}`$ as functions of $`t/t_\xi `$ for several values of $`\alpha `$ in $`d=2`$ and $`d=3`$. We obtained the best data collapse for $`\alpha =5.0\pm 0.5`$ in $`d=3`$ and $`\alpha =7\pm 1`$ in $`d=2`$, which are shown in Fig. 3a and 3b, respectively. The excellent data collapse strongly supports the above scaling assumptions.

We would like to thank Dmitry Malykhanov for the assistance with the simulations. Financial support from the German-Israeli Foundation (GIF), the Minerva Center for Mesoscopics, Fractals, and Neural Networks and the Deutsche Forschungsgemeinschaft is gratefully acknowledged.