format: text

text:
# Spectral-weight function of the ionic Hubbard model

## I  Introduction

In quasi-one dimensional materials like halogen-bridged transition metal chain complexes, conjugated polymers, organic charge transfer salts, or inorganic blue bronzes, the itinerancy of the electrons strongly competes with electron-electron and electron-phonon interactions which tend to localize the charge carriers by establishing commensurate spin or charge-density wave ground states (GS). At half filling band (BI) or Mott (MI) insulating phases are favored over the metallic state. Quantum phase transitions between the insulating phases are possible and the character of the electron excitation spectra reflect the properties of the different GS. A controversial issue is the nature of the BI-MI transition as well as whether or not only one critical point separates both insulating phases in purely electronic model Hamiltoniansresta ; fabrizio . Despite their importance in low dimensional materialsjeckel ; fehske  phonon dynamical effects will be ignored in this paper which will be concerned with the so-called one dimensional (1D) ionic Hubbard model (IHM), just a Hubbard model with a staggered potential which displaces by $`\mathrm{\Delta }`$ the energy levels of even and odd sites. It is described by the Hamiltonian

$$H=t\underset{<ij>s}{}c_{is}^+c_{js}+\frac{\mathrm{\Delta }}{2}\underset{is}{}\left(1\right)^in_{is}+U\underset{i}{}\left(n_i\frac{1}{2}\right)\left(n_i\frac{1}{2}\right),$$
(1)
where electrons of spin $`s=\pm \frac{1}{2}`$ hop with amplitude t between nearest-neighbor inequivalent sites of energy level $`\pm \frac{1}{2}\mathrm{\Delta }`$ ( $`\mathrm{\Delta }=1`$ in what follows). Electrons of opposite spin on the same site feel a repulsion $`U`$.

This model Hamiltonian, originally proposed by Nagaosa naga  and later by Egamiega  as a model ferroelectric, is ideal for studying the issue of quantum phase transitions in electron systems. On general grounds one expects a transition to occur from an ionic band insulator to a strongly-correlated Mott insulator as $`U`$ increases. Evidence for such a transition was found in exact-diagonalization calculations by Resta et al resta  predicting a metallic point at a critical value $`U=U_c`$ separating both insulating phases. This transition is signaled by a polarization jump with a sign change associated with a site parity change of the GS. The same conclusion was supported by Gidopoulos et algido , who showed the reversal of site parity above to be of magnetic origin, thus leading to a vanishing of the spin gap, $`\mathrm{\Delta }_s=0`$. Likewise, Brune et albrune  by means of the bosonization technique along with Lanczos and density matrix renormalization group (DMRG) calculations lent support to this picture.

On the other hand, Fabrizio et alfabrizio  shed new interest into the IHM in a field-theoretical bosonization analysis where they propose a new scenario for the intermediate region, $`U\mathrm{\Delta }`$, between the BI and the MI. The model should instead exhibit two quantum phase transitions: one from a BI state to a long-range bond-ordered (BO) state predicted to be in the Ising universality class and a second one from the BO to the MI state predicted to be a Kosterlitz-Thouless transition. Such transitions to BO states have recently been found in 1D Hubbard models with extended interactions (U-V) by Nakamuranaka . The BO state is a broken symmetry state in which the system becomes ferroelectric due strictly to electron-electron interactions even if all the atoms are at centers of inversion. This is also called the spontaneously dimerized insulating (SDI) phase, which, as argued by Fabrizio et alfabrizio , implies a finite spin gap, $`\mathrm{\Delta }_s`$ $`0`$. Such a SDI phase has been observed numerically by Wilkens and Martinwilkens  in a variational quantum Monte Carlo study. More precisely, they found one transition from the BI to the correlated SDI phase, but no second transition to the MI phase which is approached asymptotically by the SDI phase as $`U\mathrm{}`$. Recent DMRG calculations yasu ; qin  lead to mutually conflicting results. Thus whereas Takada and Kidoyasu  support the Resta et alresta  picture of a single first-order transition, Qin et alqin  support the Fabrizio et al fabrizio  scenario of two second-order quantum phase transitions.

In this paper we present the one-electron spectral-weight function A(k$`\omega `$) of the IHM. The calculation is done by means of a Quantum Monte Carlo (QMC) simulation in the grand canonical ensemble suplemented by an approximate effective action. Its eigenvalues and eigenvectors alow us to construct explicitly dynamic correlation functions in real time or frequency. This procedure leads to spectral functions with a reasonably high resolution. The calculation is performed at a relatively large inverse temperature, $`\beta t=10`$ ($`t=1`$ in what follows), which brings the system very close to its GS properties. Our calculated A(k$`\omega `$)clearly shows metallic character (non-vanishing density of states in the gap) at $`ka=0.5\pi `$ ($`a=1`$), between two critical U values instead of a single metallic point at the first threshold. Thus a double BI-M-MI smooth transition is predicted, the intermediate metallic regime at U$`2\mathrm{\Delta }`$ being both antiferromagnetic (AF) and dimerized

## II  Effective Hamiltonian and correlation functions

In our calculations, we use the QMC approach in the grand canonical ensemble, which has been explained in detail by Hirschhirsch ; hirsch2  and White et alwhite . As is well-known, the partition function is factorized into $`L`$ time slices of extent $`\mathrm{\Delta }\tau =\beta /L`$. For Hubbard models, the interaction $`U_in_in_i`$ is transformed into a discrete Ising fieldhirsch , $`\sigma _{il}=\pm 1`$, which depends on lattice site $`i`$ and imaginary time $`\tau =l\mathrm{\Delta }\tau `$ $`(l=1,\mathrm{}..L)`$. The non-interacting Hamiltonian for a given Ising configuration (a definite allocation of $`+1`$ and $`1`$ values for all the $`\sigma _{il}`$ is) readshirsch2

$$H(\tau )=\underset{<ij>s}{}c_{is}^+h_{ij}^s\left\{\sigma _i\left(\tau \right)\right\}c_{js}H_l$$
(2)
where

$$h_{ij}^s\left\{\sigma _i\left(\tau \right)\right\}=t_{ij}+\left(\mathrm{\Delta }\left(1\right)^i+\frac{U}{2}\mu +\lambda \alpha \sigma _{il}\right)\delta _{ij}$$
(3)
The first term on the rhs of Eq(3) is just the hopping term of the IHM, Eq(1). In the second term, $`\mu `$ is the chemical potential, $`\alpha =\pm 1`$ for up/down electron spin, and $`\lambda `$ is a non- linear function of $`U`$ given by $`\mathrm{cosh}\lambda \mathrm{\Delta }\tau =\mathrm{exp}\left(U\mathrm{\Delta }\tau /2\right)`$.In the limit $`\mathrm{\Delta }\tau 0`$ the partition function can be written as

$$Z=\underset{<\sigma >}{}Tr\text{ }\mathrm{\Pi }_l\text{ }e^{\mathrm{\Delta }\tau H_l}=\underset{<\sigma >}{}Tr\text{ }e^{\mathrm{\Delta }\tau _lH_l}=\underset{<\sigma >}{}Tr\text{ }e^{\beta (K+V)}$$
(4)
where $`_{<\sigma >}`$ means summing over the Ising configurations and $`Tr`$denotes taking the trace over the electron degrees of freedom. $`K`$ is the $`\tau `$-independent part of $`H(\tau )`$ and

$$V=\frac{\lambda }{L}\underset{il}{}\sigma _{il}\left(n_in_i\right)=\lambda \underset{i}{}<\sigma _i>\left(n_in_i\right)$$
(5)
$`<\sigma _i>=\left(1/L\right)_i\sigma _{il}`$ being the average Ising spin at the $`i`$-th site. Therefore, $`H_{eff}=K+V`$ is $`\tau `$-independent in the limit $`\mathrm{\Delta }\tau 0`$. We adopt this effective $`\tau `$-independent Hamiltonian, for doing measurements even if the number of time slices used in the calculation is finite. To minimize as much as possible size effects due to a finite $`L`$, we shall increase $`L`$ and seek for convergence. Needless to say, the update of the Ising field is done following the standard procedure of, say, Ref 15.

The corresponding Hamiltonian matrix for each spin $`s`$ (Eq(3) with $`<\sigma >`$instead of $`<\sigma _{il}>`$) is now diagonalized and, in terms of its eigenvalues $`\epsilon _{\mu s}`$ and eigenvectors $`\mu s>`$, we can construct correlation functions in real time or frecuency. For instance

$$G_{ij}^s(\omega )=\underset{\mu }{}<i\mu s>\left\{\frac{1f_{\mu s}}{\omega \epsilon _{\mu s}+i\eta }+\frac{f_{\mu s}}{\omega \epsilon _{\mu s}i\eta }\right\}<\mu sj>=G_{ij}^{s>}(\omega )+G_{ij}^{s<}(\omega )$$
(6)
for the causal Green’s functions and likewise for other quantities. As customary, $`f_{\mu s}=\left(e^{\beta \epsilon _{\mu s}}+1\right)^1`$is the Fermi-Dirac distribution and $`\eta `$ a vanishingly small imaginary part. All these quantities must now be averaged over the accepted Ising configurations to obtain the final result. It is just in this averaging process where the correlation effects are restored. Since the final result is translationally invariant, one can Fourier transform and find the spectral-weight function (SWF)

$$A_{ks}(\omega )=\frac{1}{\pi }Im\left\{G_{ks}^>(\omega )G_{ks}^<(\omega )\right\},$$
(7)
and the DOS

$$N_s(\omega )=\frac{1}{N}\underset{k^{}}{}A_{ks}(\omega )$$
(8)
which can also be obtained from (6) as $`(1/N)_iN_i(\omega )`$, $`N_i(\omega )`$ being the corresponding imaginary part of $`G_{ii}^s(\omega )`$. For our case of the IHM, this last route allows to calculate the partial DOS for even/odd atoms by simply restricting the sum over even/odd sites. As noticed earlier, this building procedure leads to a good resolution of spectral features.

## III Double insulator-metal-insulator transition

We characterize the system by its one-electron SWF, $`A_k(\omega )`$, along with the charge and spin structure factors, $`c_k`$ and $`s_k`$, defined as the Fourier transform of the static charge and spin, respectively, correlation functions. These are given as usual by $`c_{ij}=<q_iq_j>`$ and $`s_{ij}=<s_{iz}s_{jz}>`$, $`q_i=n_i+n_i`$ and $`s_{iz}=(1/2)(n_in_i)`$ being the charge and spin at the $`ith`$ site and $`\tau =0`$. This type of characterization is useful and can be of help in understanding the electronic properties of a system. For instance, the SWF tells us the presence or absence of gaps or pseudogaps as $`U`$ increases, which is of fundamental importance in establishing the insulating or metallic character and possible transitions between them of, say, the IHM. Likewise, the evolution of $`c_k`$ and $`s_k`$ with $`U`$, especially theirs peak at $`k=\pi `$, gives information about the charge or magnetic origin of the transitions.

### III.1 The spectral-weight function.

As hinted before, it is just the structure of $`A_k(\omega )`$ at $`k=\pi /2`$ which accounts for the insulating or metallic character of the IHM. Hence, Fig 1 shows this quantity for three values of $`U`$ (recall that $`t=1`$) so selected to exhibit the three fundamental regimes of the model. The calculation has been done for a chain of eighty sites with periodic boundary conditions at half filling and $`\beta =10`$. Typically 150 time slices, two hundred warm sweeps and one thousand measurements were taken. At $`U=1`$, a two-peak structure with a clear gap in between signals an insulating system. At $`U=2.25`$, both peaks have overlapped developing a third peak just at the midpoint of the gap (which becomes the Fermi level). This three-peak structure clearly shows a metallic system. Finally, at $`U=4`$ the central peak splits into two peaks leaving a gap in between. This four-peak structure signals again an insulating system which must be somehow of a different kind from the initial one (see below). As $`U`$ increases (not shown) the central gap widens and the four peaks go over into a two-peak structure as $`U\mathrm{}`$, the system becoming then a Mott insulator.

In order to estimate the approximate extension of the metallic region, Fig 2 shows the evolution of the spectral density at $`k=\pi /2`$ and $`\omega =0`$ (the midgap) as $`U`$ increases. We see a maximum at $`U_m=2.25`$ an onset threshold at $`U_11.75`$ and a falloff with an approximate endpoint at $`U_2=`$ $`2.75`$ eV leaving a zone of metallicity $`U_2U_1=1`$. Since this quantity may be size-dependent, we proceed to a finite-size scaling in order to approach the thermodynamic limit. Calculations made on chains of 8,12,16,80,120, and 160 sites show that $`U_2U_1=1`$ is roughly independent of size, but $`U_m`$ tends to go up, extrapolating to 2.50 when $`1/N0`$. Notice that only finite chains of $`N=4m`$ sites should be used since only them have $`k=\pi /2`$ as an allowed wave vector. Chains of other lengths fail to detect the metallicity. In the thermodynamic limit, of course, all chains have a wawe vector infinitessimally close to $`\pi /2`$

The three-peak structure responsible of the metallic phase disappears as soon as one departs from $`k=\pi /2`$, showing a single mainly occupied peak at $`k=(\pi /2)\delta `$ and a mainly empty peak at $`k=(\pi /2)+\delta `$ ( $`\delta =2\pi /N`$, $`N`$ being the number of sites). Hence the metallic peak should be visible in angle-resolved photoelectron spectroscopy and its inverse, but not in integrated measurements which give the DOS, Eq (8). Since the weight of that structure falls with $`N`$ (Fig 3), only a pseudogap with a sizeable DOS at the Fermi levels should be seen.

### III.2 The charge and spin structure factors

Fig 4 finally displays the evolution of the peak of $`c_k`$ and $`s_k`$ at $`k=\pi `$ for increasing U. The charge factor consistently drops from a high value at $`U=1`$, which signals a strong CDW in correspondence with a large ionicity (the fractional ionic character of the bond is $`FIC=0.32`$) which also falls down monotonically. The system, therefore starts as an ionic band insulator at small U. The spin structure factor starts close to zero, corresponding to a paramagnetic insulator, and monotonically increases with U, signaling an AF SDW of increasing amplitude. For large U, the systems evolves into a Mott insulator with a strong AF SDW and a vanishingly small CDW, i.e., zero ionicity (neutral).

In the intermediate regime, $`U_1<U_2`$, the system is both increasingly AF and neutral (decreasing ionicity, i.e., decreasing CDW). On the other hand, Table I displays a small sector of the density matrix ($`n_{ij}^s=<c_{is}^+c_{js}>`$), for $`U=2.25`$ just the top and left $`4\times 4`$ matrix, which is sufficient to show that the metallic state of this intermediate regime is dimerized. Thus $`n_{12}^s>n_{23}^s`$, etc, indicating that the sites ($`12`$), ($`34`$), and so on, form dimers. These dimers weaken for $`U>U_2`$ and tends to disappear as $`U\mathrm{}`$

## IV CONCLUDING REMARKS

Our QMC study shows three coupling regimes in the ionic Hubbard model: ($`i`$) An ionic band insulator below $`U=U_1`$.($`ii`$) An increasingly antiferromagnetic and dimerized metal which is gradually becoming neutral in the intermediate coupling regime $`U_1<U<U_2`$. And ($`iii`$) an increasingly AF and neutral insulator whose dimerization, measured by $`n_{12}n_{23}`$, is gradually dropping away for $`U>U_2`$. This last regime tends to a Mott insulator a $`U\mathrm{}`$.

A remark is here in order. Although the spectral-weight function at $`k=\pi /2`$ seems to hint a critical separation among the three coupling regimes (a threshold and an endpoint), the smooth behavior of all the integrated quantities (DOS, ionicity, energy, charge and spin structure factors, etc.) rather suggest a smooth interpolation between the three regimes which should be distinguished from true quantum phases with sharp, critical, separation among them.

To conclude, let us finally say that the results of our QMC calculations would be quite similar to those of Resta et al resta , Wilkens and Martinwilkens ,etc, if the metallic zone reduced to a single point. On the other hand, they would also be quite close to those of Fabrizio et al fabrizio  should the intermediate regime reported in this letter turn out to be insulating.

We acknowledge the financial support of the Spanish DGICYT through Project NPB 98-0683

TABLE I. Upper-left $`4\times 4`$ sector of the density matrix $`n_{ij}^s=<c_{is}^+c_{js}>`$

| 0.5477 | 0.2928 | 0.0109 | -0.0555 |
| --- | --- | --- | --- |
| 0.2928 | 0.4523 | 0.2116 | -0.0133 |
| 0.0109 | 0.2116 | 0.5477 | 0.2928 |
| 0.0555 | -0.0133 | 0.2928 | 0.4523 |


text:
# Phase transition classes in triplet and quadruplet reaction diffusion models

## I Introduction

Phase transitions in genuine nonequilibrium systems have been investigated often among reaction-diffusion (RD) type of models exhibiting absorbing states . In many cases mapping to surface growth, spin systems or stochastic cellular automata can be done. The classification of universality classes of second order transitions is still one of the most important uncompleted task. One hopes that symmetries and spatial dimensions are the most significant ingredients as in equilibrium cases, however it turned out that in many cases there is a shortage of such factors to explain novel universality classes. An important example was being investigated during the past two years that emerges at phase transitions of binary production systems  (PCPD). In these systems particle production competes with pair annihilation and single particle diffusion. If the production wins steady states with finite particle density appear in (site restricted) models with hard-core repulsion, while in unrestricted (bosonic) models the density diverges. By lowering the production/annihilation rate a doublet of absorbing states without symmetries emerges. One of such states is completely empty, the other possesses a single wandering particle. In case of site restricted systems the transition to absorbing states is continuous.

Although the nature of this transition has not completely been settled numerically and by field theory yet, an other novel class appearing in triplet production systems was proposed very recently  (TCPD). This reaction-diffusion model differs from the PCPD that for new particle generation at least three particles have to meet. It is important to note that these models do not break the DP hypothesis  —- according to which in one component systems exhibiting continuous phase transitions to single absorbing state (without extra symmetry and inhomogeneity or disorder) short ranged interactions can generate DP class transition only —- because they exhibit multiple absorbing states that are not frozen, lonely particle(s) may diffuse in them.

A phenomenologically introduced Langevin equation that exhibits real, multiplicative noise was suggested  to describe the critical behavior of reaction-diffusion models of types

$$nA(n+1)A,nAjA,$$
(1)
(with $`j<n`$ number of interacting particles) in the form

$$_t\rho (x,t)=a\rho (x,t)^n\rho (x,t)^{n+1}+D^2\rho (x,t)+\zeta (x,t),$$
(2)
with noise correlations

$$<\zeta (x,t)\zeta (x^{},t^{})>=\mathrm{\Gamma }\rho ^\mu \delta ^d(xx^{})\delta (tt^{}).$$
(3)
The classification of universality classes of nonequilibrium systems by the exponent $`\mu `$ of a multiplicative noise in the Langevin equation was suggested some time ago by Grinstein et al. . However it turned out that there may not be corresponding particle systems to real multiplicative noise cases  and an imaginary part appears as well if one derives the Langevin equation of a RD system starting from the Master equation in a proper way. This observation led Howard and Täuber to investigate systems with complex noise appearing in binary production models. Unfortunately the cases with and without occupation number restriction turned out to be different in $`d=1`$, although in $`d=2`$ this difference was found to disappear at criticality and below .

By rescaling eq.(2) one can get the corresponding mean-field critical exponents

$$\beta ^{MF}=1,\nu _{}^{MF}=n/2,\nu _{||}^{MF}=n.$$
(4)
The authors of  expect that the noise exponent should be in the range

$$1\mu n,$$
(5)
hence by simple power counting the upper critical dimension should be

$$d_c=2+\frac{42\mu }{n}.$$
(6)
This implies for a triplet processes: $`4/3d_c8/3`$ and for a quadruplet ($`n=4`$) processes: $`1d_c5/2`$.

Very recently Kockelkoren and Chaté introduced stochastic cellular automata (SCA) versions of general $`nA(n+k)A`$, $`mA(ml)A`$ type of models , where multiple particle creation on a given site is suppressed by an exponentially decreasing creation probability ($`p^{N/2}`$) of the particle number. They claim that their simulation results in 1d are in agreement with the fully occupation number restriction counterparts and set up a general table of universality classes, where as the function of $`n`$ and $`m`$ only 4 classes exist, namely the directed percolation class , the parity conserving class , the PCPD and TCPD classes.

In any case the heuristic Langevin equation with real noise assumption for RD models  should be proven for $`n>1`$. Furthermore in low dimensions topological constraints may cause different critical behavior with and without occupation number restriction . Note that in case of binary production models it had not been clear at all if the $`d_c=2`$ prediction of the bosonic field theory had also been true for site restricted systems until the numerical confirmation of . In this paper I show simulation results for lattice models with restricted site occupancy in $`d=1,2`$ with the aim of locating the upper critical dimensions and checking claims the of refs. about possible new universality classes.

## II Mean-field considerations

In this section I discuss the mean-field equation that can be set up for site restricted lattice models with general microscopic processes of the form

$$nA\stackrel{\sigma }{}(n+k)A,mA\stackrel{\lambda }{}(ml)A,$$
(7)
with $`n>1`$, $`m>1`$, $`k>0`$, $`l>0`$ and $`ml0`$. Note that this formulation is different from that of eq.(2) that is valid for coarse grained, continuous bosonic description of these reaction-diffusion systems. In this case the diffusion drops out and one can neglect the noise, hence the the competition of creation (with probability $`o\sigma `$) and annihilation or coagulation (parametrized with probability $`\lambda =1\sigma `$) is left behind

$$\frac{\rho }{t}=ak\sigma \rho ^n(1\rho )^kal(1\sigma )\rho ^m,$$
(8)
where $`\rho `$ denotes the site occupancy probability and $`a`$ is a dimension dependent coordination number. Each empty site has a probability (1-$`\rho `$) in mean-field approximation, hence the need for $`k`$ empty sites at a creation brings in a $`(1\rho )^k`$ probability factor. By expanding $`(1\rho )^k`$ and keeping the lowest order contribution one can see that for site restricted lattice systems a $`\rho ^{n+1}`$-th order term appears automatically with negative coefficient that regulates eq.(8). The steady state solution can be found analytically in many cases and may result in different, continuous or discontinuous phase transitions. Here I split the discussion of the solutions to three parts: (a) $`n=m`$, (b) $`n>m`$ and (c) $`n<m`$. In the inactive phases one expects a dynamical behavior described by the $`mA\mathrm{}`$ process, for which $`\rho t^{1/(m1)}`$ is known .

### A The $`n=m`$ symmetric case

The steady state solution in this case can be obtained by solving

$$k\sigma (1\rho )^k=l(1\sigma ),$$
(9)
where the trivial ($`\rho =0`$) solution has been factored out. For the active phase one gets

$$\rho =1\left[\frac{l}{k}\frac{1\sigma }{\sigma }\right]^{1/k},$$
(10)
which vanishes at $`\sigma _c=\frac{l}{k+l}`$ with the leading order singularity

$$\rho |\sigma \sigma _c|^{\beta ^{MF}},$$
(11)
and order parameter exponent exponent $`\beta ^{MF}=1`$. At the critical point the time dependent behavior is described by

$$\frac{\rho }{t}=2ak^2\rho ^{n+1}+O(\rho ^{n+2}).$$
(12)
that gives a leading order power-law solution

$$\rho t^{1/n}$$
(13)
hence $`\alpha _{MF}=\beta ^{MF}/\nu _{||}^{MF}=1/n`$. This was obtained from bosonic, coarse grained formulation in  too.

### B The $`n>m`$ case

In this case besides the $`\rho =0`$ absorbing state solution we can get an active state if

$$k\sigma \rho ^{nm}(1\rho )^k=l(1\sigma )$$
(14)
is satisfied. Both sides are linear functions of $`\sigma `$ such that for $`\sigma 0`$ only the $`\rho =0`$ is a solution. The left hand side is a convex function of $`\rho `$ (from above) with zeros at $`\rho =0`$ and $`\rho =1`$.

Therefore by increasing $`\sigma `$ from zero the left hand side meets the right hand side at $`\sigma _c`$, $`\rho _c>0`$ (See Fig.1(a)). If this solution is stable a first order transition takes place in the system. Note that in higher order cluster mean-field solutions, where the diffusion can play a role the transition may turn into continuous one , therefore it is important to check the type of transition for $`dd_c`$. In Section III C I shall confirm the first orderedness of such transitions for two models in 2d.

### C The $`n<m`$ case

By factoring out the trivial $`\rho =0`$ solution we are faced with the general condition for a steady state

$$k\sigma (1\rho )^k=l(1\sigma )\rho ^{mn}.$$
(15)
One can easily check that in this case the critical point is at $`\sigma _c=0`$ (see Fig.1(b)) and here the density decays with $`\alpha ^{MF}=1/(m1)`$ as in case of the $`n=1`$ branching and $`m=l`$ annihilating models showed by Cardy and Täuber  (BkARW classes). However the steady state solution for $`n>1`$ gives different $`\beta `$ exponents than those of BkARW classes, namely $`\beta ^{MF}=1/(mn)`$. This imply novel kind of critical behavior in low dimensions, that should be a subject of further investigations .

## III Simulations in two dimensions

Two dimensional simulations were performed on $`L=4001000`$ linear sized lattices with periodic boundary conditions. One Monte Carlo step (MCS) — corresponding to $`dt=1/P`$ (where $`P`$ is the number of particles) — is built up from the following processes. A particle and a random number $`x(0,1)`$ are selected randomly; if $`x<D`$ a site exchange is attempted with one of the randomly selected empty nearest neighbors (nn); if $`xD`$ $`k`$ number of new particles are created with probability $`(1p)`$ at randomly selected empty nn sites provided the number of nn particles was greater than or equal $`n`$; or if $`xD`$ $`l`$ number of particles are removed with probability $`p`$ (taking into account the $`ml0`$ condition as well). The simulations were started from fully occupied lattices and the particle density decay was measured up to $`10^610^8`$ MCS.

### A The $`3A4A`$, $`3A2A`$ symmetric triplet model

First I checked the dynamic behavior in the inactive phase for $`D=0.5`$ diffusion rate. At $`p=0.9`$ one can see the appearance of the mean-field behavior $`\rho (t)t^{1/2}`$ following $`2\times 10^6`$ MCS. By decreasing $`p`$ this scaling sets in later and later times. As Fig.2 shows for $`L=1000`$ systems with $`t_{max}=10^7`$ MCS curves with $`p0.4965`$ veer up – corresponding to the active phase — while curves with $`p0.497`$ veer down – corresponding to the absorbing state. From the $`\rho (t)`$ data I determined the effective exponents (the local slopes) defined as

$$\alpha _{eff}(t)=\frac{\mathrm{ln}\left[\rho (t)/\rho (t/m)\right]}{\mathrm{ln}(m)}$$
(16)
(where I used $`m=4`$). The critical point is estimated at $`p=0.4967(2)`$ with $`\alpha =0.33(1)`$ (for local slopes see insert of Fig.2). This value agrees with the mean-field value $`\alpha ^{MF}=1/3`$.

Density decays for several $`p`$-s in the active phase ($`0.003<ϵ=|p_cp|<0.3`$) were followed on logarithmic time scales and averaging was done over $`100`$ independent runs in a time window, which exceeds the level-off time by a decade. The steady state density in the active phase at a critical phase transition is expected to scale as

$$\rho (\mathrm{},p)|pp_c|^\beta .$$
(17)
Using the local slopes method one can get a precise estimate for $`\beta `$ as well as for the corrections to scaling

$$\beta _{eff}(ϵ_i)=\frac{\mathrm{ln}\rho (\mathrm{},ϵ_i)\mathrm{ln}\rho (\mathrm{},ϵ_{i1})}{\mathrm{ln}(ϵ_i)\mathrm{ln}(ϵ_{i1})},$$
(18)
where I used the $`p_c`$ value determined before. One can see on Fig.3 that the effective exponent for $`ϵ>0.005`$ exhibits a correction to scaling (inclined line) and tends to $`lim_{ϵ0}\beta =1.0(1)`$, which agrees with the mean-field value again. By neither the $`\alpha `$, nor the $`\beta `$ exponent can one observe logarithmic corrections suggesting $`d_c<2`$.

The density decay simulations were repeated at $`D=0.2`$, where the critical point was found at $`p_c=0.4795(1)`$ with mean-field like $`\alpha `$ exponent again.

### B The $`4A5A`$, $`4A3A`$ symmetric quadruplet model

Here simulations are much slower than in case of the triplet model, hence systems with linear size $`L=400`$ could be investigated. First I checked the dynamic behavior in the inactive phase for $`D=0.5`$. At $`p=0.9`$ a mean-field type of decay $`\rho (t)t^{1/3}`$ can be observed following $`10^6`$ MCS. As one can see on Fig.4 for $`p<0.4702`$ the density decay curves veer up, while for $`p0.4705`$ they veer down. The estimated critical point is $`p_c0.4703(1)`$. The effective exponent at $`p_c`$ extrapolates to $`\alpha 0.215(5)`$. As one can see on this graph the separatrix (critical) curve exhibits a linear shape on the $`\rho (t)t^{1/4}`$ \- $`\mathrm{ln}(t)`$ scale suggesting logarithmic corrections to scaling. Similarly the effective exponents of $`\beta `$ seem to extrapolate to $`\beta 0.71(5)`$ (Fig.3) that is very far from the mean-field value $`\beta ^{MF}=1`$. To check the possibility that a logarithmic correction can result in mean-field exponents the fitting with the lowest order correction

$$\rho (\mathrm{},p)=\left[(p_cp)(a+b\mathrm{ln}(p_cp))\right]^\beta $$
(19)
has been applied for the steady state $`\rho (\mathrm{},p)`$ data. I used the non-linear fitting of the “xmgr” graphical package with a relative error in the sum of squares with at most 0.0001. This resulted in $`\beta =1.01`$ at $`p=0.471`$ ($`a=10.8`$, $`b=6.05`$) (see insert of Fig. 3). This result in agreement with the dynamical scaling conclusion may support that the upper critical dimension for quadruplet models is $`d_c=2`$. To get more solid results further, very extensive simulations would be necessary that is beyond the scope of this study. In any case clear mean-field behavior can’t be concluded.

### C $`3A5A`$, $`2A\mathrm{}`$ and $`4A5A`$, $`2A\mathrm{}`$ hybrid models

One can find two regions in the density decay behavior by varying $`p`$ in both models. For $`p<p_c`$ steady state values are reached quickly while for $`p>p_c`$ a rapid (faster than power-law) initial density decay crosses over to $`\rho t^1`$. This is in agreement with the mean-field behavior of the $`2A\mathrm{}`$ process in 1d  dominating in the inactive phase. For the $`4A5A`$, $`2A\mathrm{}`$ quadruplet production model this threshold is at $`p_c=0.119(1)`$ (see Fig.5) where an abrupt jump is observable from $`\rho (\mathrm{})=0.833`$ to $`\rho (\mathrm{})=0`$. In case of the $`3A5A`$, $`2A\mathrm{}`$ triplet production model the threshold is at $`p_c=0.220(1)`$ with a jump from $`\rho (\mathrm{})=0.45`$ to zero.

In neither cases do we see dynamical scaling at the transition. These results are in agreement with the first order transition of the mean-field approximations given in Sect.II B for $`n>m`$.

## IV Simulations in one dimension

The simulations in one dimension were carried out on $`L=10^5`$ sized systems with periodic boundary conditions. The initial states were again fully occupied lattices, and the density of particles is followed up to $`10^6`$ MCS. An elementary MCS consists of the following processes:

$`A\mathrm{}\mathrm{}A`$ with probability D,

$`mA(ml)A`$ with probability $`p(1D)`$,

$`nA(n+k)A`$ with probability $`(1p)(1D)`$,

such that the reactions were allowed on the left or right side of the selected particle strings.

### A $`3A4A`$, $`3A2A`$ and $`3A6A`$, $`3A\mathrm{}`$ symmetric triplet models

The $`3A4A`$, $`3A2A`$ site restricted model in 1d was simulated by Park et al.  for small systems up to $`10^6`$ MCS. They concluded to find a novel kind of phase transition with the order parameter exponents $`\alpha =0.32(1)`$ and $`\beta =0.78(3)`$. For the restricted bosonic version of this model large scale simulations gave $`\alpha =0.27(1)`$ and $`\beta =0.90(5)`$ . Note that since reactive and diffusive sectors arise in this model like in PCPD, diffusion dependence or corrections to scaling may hamper to see real asymptotic behavior . Here I show extended simulation results for the strictly site restricted lattice model model with $`t_{max}=10^7`$ MCS at diffusion rate $`D=0.1`$. At the critical point the $`\alpha _{eff}(t)`$ curve exhibits a straight line shape for $`t\mathrm{}`$, while in sub(super)-critical cases $`\alpha _{eff}(t)`$ curves veer down(up) respectively. As one can see on Fig.6 following a long relaxation $`p0.3032`$ curves veer up, while $`p0.3035`$ curves veer down in the $`t\mathrm{}`$ limit. From this one can estimate $`p=p_c0.3033(1)`$ with $`\alpha =0.33(1)`$ in agreement with the results of .

By analyzing super-critical, steady state densities with the local slopes method one can read-off: $`\beta _{eff}\beta 0.95(5)`$ (see Fig.7), which is higher than the results of  and .

However one should be careful and check diffusion dependence and corrections to scaling especially because these critical exponent estimates are quite close to the mean-field values ($`\alpha ^{MF}=1/3`$, $`\beta ^{MF}=1`$) and as it was shown in Sect. III $`d_c<2`$. Since the $`p=0.303`$ and $`p=0.3035`$ $`\rho (t)`$ curves show clear curvature for large times the $`0.303<p_c<0.3035`$ conditions seems to be inevitable.

I tried to fit the steady state data in the $`0.303p0.3035`$ region by the logarithmic correction form (19) and obtained $`\beta =1.07(10)`$ at $`p_c=0.3032`$ that agrees with the mean-field value and implies $`d_c=1`$.

Just considering mean-field results, according to which $`k`$ does not play a role for $`n=m`$ models one may expect that the $`3A6A`$, $`3A\mathrm{}`$ triplet creation model exhibits the same kind of transition as the $`3A4A`$, $`3A2A`$ model. Indeed Kockelkoren and Chaté’s simulations show this . However doing lattice simulations with site restrictions it turned out that the $`3A6A`$ creation was so effective that it shifted the transition to the zero production limit ($`p=1`$) where the $`3A\mathrm{}`$ process in 1d is known to decay as $`\rho (\mathrm{ln}(t)/t)^{\frac{1}{2}}`$ . Off-critical simulations showed that $`\beta =0.33(1)`$, meaning that this transition belongs to the BkARW mean-field class. On the other hand there may be other realizations of this model, where the transition reported by Kockelkoren and Chaté is accessible.

### B The $`3A5A`$, $`2A\mathrm{}`$ model

It has been established that in $`n=1`$, $`m=l=2`$ and even $`k`$ – so called even number of offspringed branching and annihilating models (BARWe) – the parity conserving (PC) class continuous phase transition emerges . This class has also been observed in models exhibiting $`Z_2`$ symmetric absorbing states, where the domain walls separating ordered phases follow BARWe dynamics . This class was originally called parity conserving, owing to the conservation law that made it different from the robust DP class. However it turned out that in $`Z_2`$ symmetric models this conservation is not enough . Furthermore in binary spreading models this conservation was found to be irrelevant . Therefore it is still an open question whether parity conservation is relevant in other models than in BARW types.

I investigated the phase transition of the triplet production $`3A5A`$, $`2A\mathrm{}`$ model (with explicit particle diffusion) possessing parity conservation. As I showed in Sect. III C in two dimensions this system exhibits a first order transition in agreement with the mean-field results. This first order mean-field behavior does not give a direct hint on the type of phase transition in 1d. Kockelkoren and Chaté’s simulations on the one dimensional, suppressed bosonic cellular automaton version of this model shows simple DP class density decay . However if we consider the space-time evolution we see very non-DP like spatio-temporal pattern (see Fig.8). This pattern resembles much more to those of the PCPD class models, where compact domains separated by clouds of lonely wandering particles occur. Of course such qualitative judgment on the universal behavior is not enough but has been found to be quite successful in case of binary production systems .

The density decay simulations at $`D=0.8`$ and $`D=0.2`$ have been analyzed by the local slopes method see Fig.9. At $`D=0.8`$ the critical point is estimated at $`p_c^H=0.4629(3)`$ and the corresponding effective exponent tends to $`\alpha ^H=0.24(1)`$. At $`D=0.2`$-s the critical point is at $`p_c^L=0.2240(3)`$, and the local exponent seems to extrapolate to $`\alpha ^L=0.28(1)`$. Such small difference between the high and low $`D`$ results can also be observed by analyzing the steady state results: $`\beta ^H=0.43(3)`$ versus $`\beta ^L=0.63(3)`$. These exponent estimates are far from the 1+1 d DP values ($`\alpha =0.159464(6)`$, $`\beta =0.276486(8)`$ ), hence the claim of Kockelkoren and Chaté for the critical behavior of $`n>m`$ models is questionable.

On the other hand the diffusion dependence of the critical exponents is a challenge and has been observed in the binary production PCPD model . In  it was shown that assuming logarithmic corrections to scaling – that is quite common in 1d models – a single universality class can be supported numerically. Therefore here again I have investigated the possibility of the collapse of the high and low $`D`$ exponents. Assuming the same kind of logarithmic correction forms as in

$$\left[[a+b\mathrm{ln}(t)]/t\right]^\alpha $$
(20)
I have found a consistent set of exponents both for $`D=0.2`$ and $`D=0.8`$ (see Table I and insert of Fig.9). For the data analysis I used non-linear fitting of the program “xmgr” package, with a relative error in the sum of squares with at most 0.001.

$$\alpha =0.22(1),\beta =0.60(1),$$
(21)
with the critical thresholds: $`p_c^H=0.4627(1)`$, $`p_c^L=0.2240(1)`$. These exponents suggest a distinct universality class from the known ones .

### C $`4A5A`$, $`4A\mathrm{}`$ and $`4A5A`$, $`2A\mathrm{}`$ quadruplet models

Two dimensional simulations (Sect.III), showed that for $`n=m=4`$ symmetric quadruplet models $`d_c=2`$. Simulations in the corresponding suppressed bosonic SCA  with $`n=4`$ and $`1m4`$ located the phase transition at zero production rate. Here I show that in the one dimensional $`4A5A`$, $`4A\mathrm{}`$ and $`4A5A`$, $`2A\mathrm{}`$ site restricted models continuous phase transitions with $`p<1`$ and with non-trivial exponents can be found. The density decay was followed up to $`t=10^6`$ MCS and the critical point was located by the local-slopes method (see Fig.10) at $`p=0.9028(1)`$ for $`D=0.3`$. The corresponding exponent can be estimated as $`\alpha =0.27(1)`$. For $`D=0.05`$ one gets $`p_c=0.9605(3)`$ with $`\alpha =0.28(1)`$, so one can not see diffusion dependence here. Analyzing off-critical data with the local slopes method (18) one gets $`\beta =0.48(2)`$ (see Fig.7).

In accordance with these results simulations for the $`4A5A`$, $`2A\mathrm{}`$ model at $`D=0.2`$ and $`D=0.8`$ diffusion rates resulted in $`p_c(0.2)=0.53185(5)`$, $`p_c(0.8)=0.5742(1)`$ with $`\alpha =0.27(1)`$ and $`\beta =0.48(2)`$ exponents (see Fig.7). As we can see critical exponent data for quadruplet models are robust and no diffusion dependence has been found. Furthermore critical space-time plots are very similar to that of the PCPD model.

## V Conclusions

In this paper I investigated the phase transitions of general $`nA(n+k)A`$, $`mA(ml)A`$ reaction type of models with explicit single particle diffusion on occupation number restricted lattices in one and two dimensions. I showed that mean-field solution for $`n=m`$ symmetric cases results in universality classes characterized by the exponents $`\alpha =1/n`$, $`\beta =1`$. I determined the upper critical dimensions for the triplet and quadruplet cases by simulations. For $`n=3`$ high precision simulations show mean-field type of criticality with logarithmic corrections meaning $`d_c=1`$. This result is in contradiction with the simulations of  and  and with the analytical form for $`d_c(n)`$ derived from a phenomenological Langevin equation. In case of my site restricted realization of the one dimensional $`3A6A`$, $`3A\mathrm{}`$ model the phase transition point is shifted to zero production rate and is continuous, BkARW mean-field type. This is in contradiction with the findings of  for an other stochastic cellular automaton realization of this model. For $`n=4`$ the upper critical dimension was located at $`d_c=2`$ opening up the possibility for non-trivial critical behavior in $`d=1`$. Indeed two versions of such quadruplet models were shown to exhibit robust, novel type of critical transition in one dimension.

For $`n>m`$ the mean-field approximation gives first order transition that was observed by simulations for two ($`n=3,4`$) models in $`d=2`$. On the other hand numerical evidence was given that the parity conserving model $`3A5A`$, $`2A\mathrm{}`$ in 1d exhibits a non-PC type of critical behavior with logarithmic corrections by varying the diffusion rate. This transition does not fit in the universality class scheme suggested by .

Finally I showed that for $`n<m`$ models the mean-field approximations result in new classes featured by $`\alpha ^{MF}=1/(m1)`$ and $`\beta ^{MF}=1/(mn)`$. Such kind of models should be subject of further studies.

The presented mean-field and simulation results show that the universal behavior of such low-dimensional reaction-diffusion models is rich and the table of universality classes given by  is not valid for 1d, fully site restricted systems. Perhaps the strict site restriction plays an important role that causes the differences. Field theoretical (possibly fermionic) treatment that starts from the master equation should be set up to determine at least the analytical form of $`d_c(n)`$ for $`n=m>2`$ models.

Acknowledgements:
I thank H. Chaté, H. Hinrichsen and U. Täuber for useful communications and N. Menyhárd for her comments to the manuscript. Support from Hungarian research funds OTKA (Grant No. T-25286), Bolyai (Grant No. BO/00142/99) and IKTA (Project No. 00111/2000) is acknowledged. The simulations were performed on the parallel cluster of SZTAKI and on the supercomputer of NIIF Hungary within the scope of the DEMOGRID project.


text:
# A simple construction of the fractional Brownian motion.

## 1. Introduction

The fractional Brownian motion appears to be a very natural object, for its three fundamental and characteristic features of being a continuous Gaussian process, self-similar and with stationary increments. By self-similar, we mean that there exists a real number $`H]0,1[`$, such that the finite-dimensional distributions of $`\{T^HB_H(Tt),t0\}`$ do not depend on $`T`$. The parameter $`H`$ is called the Hurst parameter or the index of self-similarity. Under these conditions, it is not hard to check that the covariance function must have the following form:

$$E[B_H(s)B_H(t)]=\frac{1}{2}(s^{2H}+t^{2H}|st|^{2H}).$$

A major consequence of this fact is that if we introduce the sequence of increments $`G_H(j)=B_H(j)B_H(j1),j=1,2,\mathrm{}`$, also called fractional Gaussian noise, we note that they are strongly correlated (for $`H\frac{1}{2}`$). More precisely,

$$E[G_H(j)G_H(j+k)]\underset{k\mathrm{}}{}H(2H1)k^{2H2}(1)$$

We denote the two radically different behaviours on both sides of $`\frac{1}{2}`$: for $`H<1/2`$ the increments are all negatively correlated, which corresponds to a chaotic behavior, whereas for $`H>1/2`$ the positive correlation between the increments corresponds to a more disciplined behaviour. We refer to Chapter 7 of Samorodnitsky and Taqqu  for more information about fractional Brownian motion.

There is a well-known representation introduced by Mandelbrot and Van Ness  of the fractional Brownian motion as an integral of a kernel function with respect to the usual Brownian motion. The method of approximation that uses the discretization of this integral leads to an algorithm obliging to store a lot of data in memory, and to deal with non smooth functions in the case $`0<H<\frac{1}{2}`$. This inconvenience is nicely discussed by Carmona and Coutin in , and they propose a way to reduce it.

The construction we propose is based on correlated random walks: it consists in discrete processes such that the law of each move is ruled by the value of the previous move. We refer to  and the references therein for more information about these processes. We note that the decay of correlation for such processes is exponential, but that a mixture of these processes, adapted to the value of the index of self-similarity, leads to walks whose correlations satisfy (1). The superposition of a large number of such walks yields a discrete Gaussian process whose correlations fulfill the conditions of Taqqu , so that its scaling limit is the fractional Brownian motion. We have to distinguish both cases $`\frac{1}{2}H<1`$ and $`0<H<\frac{1}{2}`$. In the first case, the equivalence relation (1) is the only condition to check. The case $`0<H<\frac{1}{2}`$ is more delicate: a compensation relation between the correlations has to be satisfied simultanously. That leads us to introduce two different types of correlated random walks for both cases.

This construction reminds the constructions using renewal processes of Mandelbrot  and Taqqu and Levy ,. It reminds even more the construction coming from traffic modeling by Taqqu, Willinger and Sherman . To our knowledge, all these methods are restricted to the case $`\frac{1}{2}H<1`$. In that case, we shall discuss the differences between our construction and . Finally, we want to mention that random media has already been used by Kesten and Spitzer  to get convergences towards certain self-similar processes with index of self-similarity bigger than $`\frac{1}{2}`$.

In the last part, we discuss the speed of the convergence stated in the previous parts. We indicate how to deduce a quite fast algorithm that does not ask to keep a lot of datas in memory.

## 2. The case $`1/2H<1`$

### 2.1. The correlated random walk

We introduce first our basic tool: the correlated random walk with persistence $`p`$. It is a process evolving on $``$ by jumps of +1 or -1, whose probability of making the same jump as the previous one is $`p`$:

###### Definition 1.

For any $`p[0,1]`$, the correlated random walk $`X^p`$ with persistence $`p`$ is a $``$-valued discrete process, such that :

_ $`X_0^p=0`$, $`P(X_1^p=1)=1/2`$, $`P(X_1^p=1)=1/2`$.

_ $`n1,\epsilon _n^p:=X_n^pX_{n1}^p`$ equals 1 or -1 a.s.

_ $`n1,P(\epsilon _{n+1}^p=\epsilon _n^p|\sigma (X_k^p,0kn))=p.`$

This process is not Markovian, but if we define a state as the position of the process on $``$, coupled with the sign of its last jump, we have to deal with a Markov process on $`\times \{1,1\}`$. In fact, this process consists in alternate falls and rises of i.i.d. geometric length with parameter $`p`$.

We can compute the correlations between two steps distant from $`n`$:

###### Proposition 1.

$`m1,n0,`$ $`E[\epsilon _m^p\epsilon _{m+n}^p]=(2p1)^n`$

Proof: $`n1`$,

$`E[\epsilon _{n+1}^p|\sigma (X_k,0kn)]=`$

$`E[\epsilon _{n+1}^p1_{\epsilon _n^p=1}|\sigma (X_k,0kn)]+E[\epsilon _{n+1}^p1_{\epsilon _n^p=1}|\sigma (X_k,0kn)]=`$

$`(2p1)1_{\epsilon _n^p=1}(2p1)1_{\epsilon _n^p=1}=(2p1)\epsilon _n^p`$.

Consequently, conditioning by $`\sigma (X_k,0km+n)`$, we get

$`m1,n0,`$ $`E[\epsilon _{m+n+1}^p\epsilon _m^p]=(2p1)E[\epsilon _{m+n}^p\epsilon _m^p]`$.

The result is then obtained by recurrence. ∎

We now introduce an extra randomness in the persistence. We first denote by $`P^p`$ the law of $`X^p`$ for a given $`p`$. Now, considering a probability measure $`\mu `$ on $`[0,1]`$, we will call $`P^\mu `$, the annealed law of the correlated walk associated to $`\mu `$, i.e. the measure on $`^{}`$ defined by $`P^\mu :=_0^1P^p𝑑\mu (p)`$.

Remark: unlike the situation in , the persistence does not depend on the level. Only one coin toss, according to $`\mu `$, decides for the whole environment.

Let $`X^\mu `$ be a process of law $`P^\mu `$. Let us now introduce the notation $`\epsilon _n^\mu :=X_n^\mu X_{n1}^\mu `$. From Proposition 1 we get the straightforward result:

###### Proposition 2.

$`m1,n0,`$ $`E[\epsilon _m^\mu \epsilon _{m+n}^\mu ]=_0^1(2p1)^n𝑑\mu (p).`$

### 2.2. Statement and proof of the result

The goal now is to introduce a probability measure $`\mu `$ leading to the same equivalent as (1), mentioned in the introduction, so that by taking the average over a large number of trajectories, we approximate a discrete Gaussian process having the same properties as in , whose scaling limit is the fractional Brownian motion:

###### Theorem 1.

Let $`H]1/2,1[`$.

Denote by $`\mu ^H`$ the probability on $`[\frac{1}{2},1]`$ with density $`(1H)2^{32H}(1p)^{12H}`$.

Let $`(X^{\mu ^H,i})_{i1}`$ be a sequence of independent processes of law $`P^{\mu ^H}`$,

$$^𝒟\underset{N\mathrm{}}{lim}\underset{M\mathrm{}}{lim}c_H\frac{X_{[Nt]}^{\mu ^H,1}+\mathrm{}+X_{[Nt]}^{\mu ^H,M}}{N^H\sqrt{M}}=B_H(t)$$

with $`c_H=\sqrt{\frac{H(2H1)}{\mathrm{\Gamma }(32H)}}`$,

$``$ means convergence in the sense of the finite-dimensional distributions, and $`^𝒟`$ means convergence in the sense of the weak convergence in the Skorohod topology on $`D[0,1]`$, the space of cadlag functions on .

Proof: The central limit theorem implies that $`lim_M\mathrm{}\frac{X_k^{\mu ^H,1}+\mathrm{}+X_k^{\mu ^H,M}}{\sqrt{M}}`$ is a discrete centered Gaussian process $`(Y_k^H)_{k1}`$, with stationary increments $`G_k^H:=Y_{k+1}^HY_k^H`$ with $`E[G_k^H]=0`$, $`E[(G_k^H)^2]=1`$ and

$$i,n0,r(n):=E[G_i^HG_{i+n}^H]=(22H)2^{22H}_{\frac{1}{2}}^1(2u1)^n(1u)^{12H}𝑑u$$

$`\begin{array}{cc}\hfill r(n)& =(22H){\displaystyle _0^1}v^n(1v)^{12H}𝑑v\hfill \\ & =(22H){\displaystyle \frac{\mathrm{\Gamma }(n+1)\mathrm{\Gamma }(22H)}{\mathrm{\Gamma }(n+32H)}}+O({\displaystyle \frac{1}{n}})\hfill \\ & \underset{n\mathrm{}}{}\mathrm{\Gamma }(32H){\displaystyle \frac{1}{n^{22H}}}={\displaystyle \frac{1}{c_H^2}}{\displaystyle \frac{H(2H1)}{n^{22H}}}\hfill \end{array}`$

So that,

$`\begin{array}{cc}\hfill E[c_H^2(G_1^H+\mathrm{}+G_N^H)^2]& =c_H^2{\displaystyle \underset{i=1}{\overset{N}{}}}{\displaystyle \underset{j=1}{\overset{N}{}}}r(|ij|)\hfill \\ & =c_H^2(r(0)+{\displaystyle \underset{i=1}{\overset{N1}{}}}[r(0)+2{\displaystyle \underset{k=1}{\overset{i}{}}}r(k)])\hfill \\ & \underset{n\mathrm{}}{}N^{2H}\hfill \end{array}`$

(the last step consists simply in two successive comparisons between sums and integrals). A direct application of  (lemma 5.1) allows to conclude. ∎

We can give also an analog statement for $`H=1/2`$:

###### Theorem 2.

Denote by $`\mu ^{\frac{1}{2}}`$ the uniform probability on $`[\frac{1}{2},1]`$.

Let $`(X^{\mu ^{\frac{1}{2}},i})_{i1}`$ be a sequence of independent processes of law $`P^{\mu ^{\frac{1}{2}}}`$,

$$^𝒟\underset{N\mathrm{}}{lim}\underset{M\mathrm{}}{lim}c_{\frac{1}{2}}\frac{X_{[Nt]}^{\mu ^{\frac{1}{2}},1}+\mathrm{}+X_{[Nt]}^{\mu ^{\frac{1}{2}},M}}{\sqrt{N\mathrm{log}N}\sqrt{M}}=B(t)$$

where $`B`$ is the classical Brownian motion, and $`c_{\frac{1}{2}}=\frac{1}{\sqrt{2}}`$.

Proof: The scheme is the same as in Theorem 1. The difference here is that $`r(n)=2_{\frac{1}{2}}^1(2u1)^n𝑑u=\frac{1}{(n+1)}`$.

So that, $`r(0)+_{i=1}^{N1}[r(0)+2_{k=1}^ir(k)]\underset{n\mathrm{}}{}2N\mathrm{log}N`$.

We conclude again, applying  (lemma 5.1). ∎

Remark: The order of the limits in both theorems is of big importance: the limit in the reverse order would bring 0, as far as for any fixed $`p`$, a correlated random walk satisfies a central limit theorem with normalization $`\sqrt{N}=o(\sqrt{N\mathrm{log}N})`$ (see ).

We want now to compare previous theorems with the result of : in , the limit theorem consists also in taking the scaling limit of the average over a large number (tending to infinity) of i.i.d. copies of processes which are a succession of falls and rises. In the case of , falls and rises are all independent with infinite variance.

In our setting, even if the lengths of falls and rises have finite variance under each $`P^p`$ (they are geometric), the laws under $`P^{\mu ^H}`$ of the falls and the rises (which are all the same) have infinite variance: indeed, if we denote by $`L`$ the length of a rise,

$$\begin{array}{cc}\hfill P^{\mu ^H}(Ln)& =_0^1P^p(Ln)𝑑\mu ^H(p)=_0^1p^n𝑑\mu ^H(p)\hfill \\ & =(22H)_0^1p^n(1p)^{12H}𝑑p\underset{n\mathrm{}}{}\mathrm{\Gamma }(32H)\frac{1}{n^{22H}}\hfill \end{array}$$

(for $`H=1/2`$, we get $`P^{\mu ^{\frac{1}{2}}}(Ln)\underset{n\mathrm{}}{}{\displaystyle \frac{1}{n}}`$)

and we get the same kind of tail as in .

The difference lies in the fact that the falls and the rises are not independent: indeed, if the first rise is short, it probably means that the environment $`p`$ is small, so that the following fall will be probably short also. More precisely, if $`L_1`$ and $`L_2`$ denote the lengths of the first rise (resp. fall) and of its following fall (resp. rise),

$$\begin{array}{cc}\hfill P^{\mu ^H}(L_1n,L_2m)& =E_{\mu ^H}[P^p(L_1n,L_2m)]=E_{\mu ^H}[p^{n+m}]\hfill \\ & \underset{n\mathrm{}}{}\mathrm{\Gamma }(32H)\frac{1}{(n+m)^{22H}}\hfill \end{array}$$

$`L_1`$ and $`L_2`$ are therefore not independent.

Note: Theorems 1 and 2 can be extended to any probability measure with moments equivalent to $`\frac{1}{n^{22H}L(n)}`$, where $`L`$ is a slowly varying function. The additional arguments can be found in  and rely mainly on Karamata’s theorem in order to replace the naive comparison between sums and integrals we used at the end of the proof. We chose $`\mu ^H`$ in our statement for the simple formula it yields for $`c_H`$ and because it is the law of $`1\frac{U^{\frac{1}{22H}}}{2}`$ where $`U`$ is uniform on $`[0,1]`$, which makes it easy to simulate. At the end of the article, we will discuss the practical interest of other measures.

## 3. The case $`0<H<1/2`$

### 3.1. The alternating correlated random walk

We first remark that the correlated random walks of section 2, cannot provide negative correlations for the increments, at least for increments separated by an even time interval. The best we can hope is to get an alternate sign for the correlations. In order to get a process with always negative correlations (except for variances), we will consider the sequence of the sum of two consecutive increments. More precisely, if we consider, with the notations of section 2, the sequence $`(\epsilon _{2n+1}^p+\epsilon _{2n+2}^p)_{n0}`$, for any $`p`$ less than $`1/2`$, we get indeed a sequence of negatively correlated variables and it is also possible to exhibit a probability on $`[0,1]`$ such that the equivalence relation (1) mentioned in the introduction will be satisfied.

But, in the case $`0<H<1/2`$, this condition alone does not ensure a scaling limit for the scale $`N^H`$ (which is that time smaller than $`\sqrt{N}`$). It has to be allied to a compensation relation between all the correlations. We refer to  (section 5) for the statement of this condition, and we will explicit it further.

It is the reason why we have to introduce a different kind of walk, we will call alternating correlated random walk with persistence $`p`$. It is a process evolving on $``$ by jumps of +1 or -1, whose probability of making the same jump as the previous one is alternately $`p`$ and 0. In other words, one jump over two is the opposite of the previous one:

###### Definition 2.

For any $`p[0,1]`$, the alternating correlated random walk $`\stackrel{~}{X}^p`$ with persistence $`p`$, is a $``$-valued discrete process, such that :

_ $`\stackrel{~}{X}_0^p=0`$, $`P(\stackrel{~}{X}_1^p=1)=1/2`$, $`P(\stackrel{~}{X}_1^p=1)=1/2`$.

_ $`n1,\stackrel{~}{\epsilon }_n^p:=\stackrel{~}{X}_n^p\stackrel{~}{X}_{n1}^p`$ equals 1 or -1 a.s.

_ $`n1,P(\stackrel{~}{\epsilon }_{2n}^p=\stackrel{~}{\epsilon }_{2n1}^p|\sigma (\stackrel{~}{X}_k^p,0k2n1))=p.`$

_ $`n1,\stackrel{~}{\epsilon }_{2n+1}^p=\stackrel{~}{\epsilon }_{2n}^p`$.

As suggested in the introduction of this section, we will be actually interested by the process

$$Y_n^p:=\frac{\stackrel{~}{X}_{2n}^p}{2\sqrt{p}}$$

(the importance of this normalization will appear later). The trajectories of this process take only two values, which are either $`1/\sqrt{p}`$ and 0 or 0 and $`1/\sqrt{p}`$. The successive lengths of the time intervals during which the process stays on each value are independent and geometric with parameter $`1p`$.

We compute now the correlations of the increments of $`Y_n^p`$ i.e. of $`\delta _n^p:=Y_n^pY_{n1}^p=\frac{1}{2\sqrt{p}}(\stackrel{~}{\epsilon }_{2n1}^p+\stackrel{~}{\epsilon }_{2n}^p)`$ for $`n1`$.

###### Proposition 3.

_ $`m1,`$ $`E[(\delta _m^p)^2]=1`$

_ $`m1,n1,`$ $`E[\delta _m^p\delta _{m+n}^p]=p(12p)^{n1}`$

Proof: As in Proposition 1, everything is based on the following facts:

$`n1`$,

_ $`E[\stackrel{~}{\epsilon }_{2n}^p|\sigma (\stackrel{~}{X}_k^p,0k2n1)]=(2p1)\stackrel{~}{\epsilon }_{2n1}^p`$,

_ $`E[\stackrel{~}{\epsilon }_{2n+1}^p|\sigma (\stackrel{~}{X}_k^p,0k2n)]=\stackrel{~}{\epsilon }_{2n}^p`$,

$`m1,`$

$`\begin{array}{cc}\hfill E[(\delta _m^p)^2]& =\frac{1}{4p}(E[(\stackrel{~}{\epsilon }_{2m1}^p)^2]+E[(\stackrel{~}{\epsilon }_{2m}^p)^2]+2E[\stackrel{~}{\epsilon }_{2m1}^p\stackrel{~}{\epsilon }_{2m}^p])\hfill \\ & =\frac{1}{4p}(2+2(2p1))\hfill \\ & =1\hfill \end{array}`$

$`m1,n1,`$ by successive conditionings,

$`\begin{array}{cc}\hfill E[\delta _m^p\delta _{m+n}^p]& =\frac{1}{4p}E[(\stackrel{~}{\epsilon }_{2m1}^p+\stackrel{~}{\epsilon }_{2m}^p)(\stackrel{~}{\epsilon }_{2m+2n1}^p+\stackrel{~}{\epsilon }_{2m+2n}^p)]\hfill \\ & =\frac{1}{4p}(12p)^{n1}E[(\stackrel{~}{\epsilon }_{2m1}^p+\stackrel{~}{\epsilon }_{2m}^p)(\stackrel{~}{\epsilon }_{2m+1}^p+\stackrel{~}{\epsilon }_{2m+2}^p)]\hfill \\ & =\frac{1}{4p}(12p)^{n1}(12(2p1)(2p1)^2)\hfill \\ & =p(12p)^{n1}\hfill \end{array}`$

We can already note the following (compensation) relation:

$$E[(\delta _m^p)^2]+2\underset{n1}{}E[\delta _m^p\delta _{m+n}^p]=0(2)$$

Again we introduce an extra randomness in the persistence. We first denote by $`Q^p`$ the law of $`Y^p`$ for a given $`p`$. Now, considering a probability measure $`\mu `$ on $`[0,1]`$, we will call $`Q^\mu `$, the annealed law of the correlated walk associated to $`\mu `$, i.e. the measure on $`^{}`$ defined by $`dQ^\mu :=_0^1Q^p𝑑\mu (p)`$.

Let $`Y^\mu `$ be a process of law $`Q^\mu `$. We introduce $`\delta _n^\mu :=Y_n^\mu Y_{n1}^\mu `$. From Proposition 3 we get the straightforward result:

###### Proposition 4.

_ $`m1,`$ $`E[(\delta _m^\mu )^2]=1.`$

_ $`m1,n1,`$ $`E[\delta _m^\mu \delta _{m+n}^\mu ]=_0^1p(12p)^{n1}𝑑\mu (p).`$

Note: As the compensation relation (2) is satisfied ”$`p`$ by $`p`$”, it remains true for the annealed correlations.

The goal now is to introduce a probability measure $`\mu `$ leading to the same equivalent as (1), mentioned in the introduction.

### 3.2. Statement and proof of the result

We proceed as in previous section:

###### Theorem 3.

Let $`H]0,1/2[`$.

Denote by $`\mu ^H`$ the probability on $`[0,\frac{1}{2}]`$ with density $`(12H)2^{12H}p^{2H}`$.

Let $`(Y^{\mu ^H,i})_{i1}`$ be a sequence of independent processes of law $`P^{\mu ^H}`$,

$$^𝒟\underset{N\mathrm{}}{lim}\underset{M\mathrm{}}{lim}c_H\frac{Y_{[Nt]}^{\mu ^H,1}+\mathrm{}+Y_{[Nt]}^{\mu ^H,M}}{N^H\sqrt{M}}=B_H(t)$$

with $`c_H=\sqrt{\frac{2H}{\mathrm{\Gamma }(22H)}}`$.

Proof: The central limit theorem implies that $`lim_M\mathrm{}\frac{Y_k^{\mu ^H,1}+\mathrm{}+Y_k^{\mu ^H,M}}{\sqrt{M}}`$ is a discrete centered Gaussian process $`(Z_k^H)_{k1}`$, with stationary increments $`W_k^H:=Z_{k+1}^HZ_k^H`$ with $`E[W_k^H]=0`$, $`E[(W_k^H)^2]=1`$ and

$$i,n0,r(n):=E[W_i^HW_{i+n}^H]=(12H)2^{12H}_0^{\frac{1}{2}}(12u)^nu^{12H}𝑑u$$

$`\begin{array}{cc}\hfill r(n)& ={\displaystyle \frac{(12H)}{2}}{\displaystyle _0^1}(1v)^nv^{12H}𝑑v\hfill \\ & ={\displaystyle \frac{(12H)}{2}}{\displaystyle \frac{\mathrm{\Gamma }(n+1)\mathrm{\Gamma }(22H)}{\mathrm{\Gamma }(n+32H)}}\hfill \\ & \underset{n\mathrm{}}{}{\displaystyle \frac{(12H)}{2}}\mathrm{\Gamma }(22H){\displaystyle \frac{1}{n^{22H}}}={\displaystyle \frac{1}{c_H^2}}{\displaystyle \frac{H(2H1)}{n^{22H}}}\hfill \end{array}`$

So that,

$`\begin{array}{cc}\hfill E[c_H^2(G_1^H+\mathrm{}+G_N^H)^2]& =c_H^2{\displaystyle \underset{i=1}{\overset{N}{}}}{\displaystyle \underset{j=1}{\overset{N}{}}}r(|ij|)\hfill \\ & =c_H^2(r(0)+{\displaystyle \underset{i=1}{\overset{N1}{}}}[r(0)+2{\displaystyle \underset{k=1}{\overset{i}{}}}r(k)])\hfill \\ & =c_H^2(r(0)2{\displaystyle \underset{i=1}{\overset{N1}{}}}{\displaystyle \underset{k=i+1}{\overset{\mathrm{}}{}}}r(k))\hfill \\ & \underset{n\mathrm{}}{}N^{2H}\hfill \end{array}`$

The last equality comes from the compensation relation (2), and the last step consists simply in two successive comparisons between sums and integrals. A direct application of  (lemma 5.1) allows to conclude. ∎

The process $`Y_n^p`$ reminds strikingly the On/Off processes of the traffic modeling theory, described in , or the $`\mathrm{"}V\mathrm{"}`$-renewal process of  which would have alternate rewards. But our attempt to make this model work here, failed probably because it is not clear how to translate to this model the normalization of $`\stackrel{~}{X}_{2n}^p`$ by $`\sqrt{p}`$, which is crucial in our construction because its absence would yield an infinite measure for $`\mu ^H`$.

Remark 1: As in previous section, $`\mu ^H`$ is not the only convenient measure. We remark that it is the law of $`\frac{U^{\frac{1}{12H}}}{2}`$, where $`U`$ is uniform on $`[0,1]`$.

Remark 2: Other normalizations of $`\stackrel{~}{X}_{2n}^p`$ by powers of $`p`$ can be considered, but then the measure $`\mu ^H`$ has to be multiplied by a power of $`p`$ also.

## 4. Practical aspects

As mentioned above, taking the limit in reverse order, yields trivial processes. The question is: for a given number of steps $`N`$, what kind of $`M`$ are we supposed to take in order to approximate the right limit ? We restrict our study to the law of $`B_H(1)`$. We base our study on Berry-Esseen’s inequality, applied to the sequence of i.i.d. variables with law $`c_HX_N^{\mu ^H}/N^H`$ (resp. $`c_HY_N^{\mu ^H}/N^H`$) for $`H1/2`$ (resp. $`H<1/2`$).

What remains to compute, is an upper bound for the third moment of the absolute value of these variables.

We begin with the case $`H>\frac{1}{2}`$.

###### Proposition 5.

For $`H>1/2`$, for $`N`$ large enough,

$$E[(\frac{c_H|X_N^{\mu ^H}|}{N^H})^3]D_HN^{1H},$$

with $`D_H=\sqrt{\frac{6(2H1)}{(H+1)(2H+1)}}\times c_H.`$

Nota Bene: we express $`D_H`$ in terms of $`c_H`$, in order to have a formula that works for other measures, as it will be useful in the following.

Proof: We omit here the superscipts $`\mu ^H`$ in the variables $`X_n^{\mu ^H}`$’s and $`\epsilon _n^{\mu ^H}`$’s. $`X_N=_{k=1}^N\epsilon _k`$, with:

_ $`\epsilon _n`$ are Bernoulli(1/2),

_ $`c_H^2Cov(\epsilon _k,\epsilon _l)=r(|kl|)`$, with $`r(n)_{\mathrm{}}\frac{H(2H1)}{n^{22H}}`$.

We use Cauchy-Schwarz inequality to get:

$$\begin{array}{cc}\hfill E[(\frac{c_H|X_N|}{N^H})^3]& E[(\frac{c_H|X_N|}{N^H})^2]^{\frac{1}{2}}E[(\frac{c_H|X_N|}{N^H})^4]^{\frac{1}{2}}\hfill \\ & \frac{c_H^2}{N^{2H}}(\underset{1i_kN}{}_{1k4}E[\epsilon _{i_1}\epsilon _{i_2}\epsilon _{i_3}\epsilon _{i_4}])^{1/2}\hfill \end{array}$$

(using that the variance converges to 1 for large $`N`$).

Assume $`i_4i_3i_2i_1`$, we get, as in Proposition 1, by successive conditionings:

$$E[\epsilon _{i_1}^p\epsilon _{i_2}^p\epsilon _{i_3}^p\epsilon _{i_4}^p]=(2p1)^{(i_1i_2)+(i_3i_4)}$$

So that,

$$E[\epsilon _{i_1}\epsilon _{i_2}\epsilon _{i_3}\epsilon _{i_4}]=r((i_1i_2)+(i_3i_4))$$

We see by using the equivalence between sums and integrals, as we did at the end of the proof of Theorem 1, that the sum is equivalent to:

$$4!\times _0^N_0^{x_1}_0^{x_2}_0^{x_3}r((x_1x_2)+(x_3x_4))𝑑x_4𝑑x_3𝑑x_2𝑑x_1$$

$$\frac{4!}{c_H^2}\times H(2H1)_0^N_0^{x_1}_0^{x_2}_0^{x_3}((x_1x_2)+(x_3x_4))^{2H2}𝑑x_4𝑑x_3𝑑x_2𝑑x_1$$

$$=\frac{6}{c_H^2}\frac{(2H1)N^{2H+2}}{(H+1)(2H+1)}$$

Now applying Berry-Esseen’s inequality, the error of the distribution function of the marginal at time 1, is dominated by $`0.65\times D_H\frac{N^{1H}}{\sqrt{M}}`$ (using the constant 0.65 of , as far as the third moment is much bigger than the power $`\frac{3}{2}`$ of the variance).

We deduce that Theorem 1 remains true as soon as $`M(N)\mathrm{}`$ as $`N\mathrm{}`$ at a faster rate than $`N^{22H}`$:

###### Corollary 1.

Let $`M`$ be a function on the integers such that $`\frac{M(N)}{N^{22H}}`$ tends to $`\mathrm{}`$,

$$^𝒟\underset{N\mathrm{}}{lim}c_H\frac{X_{[Nt]}^{\mu ^H,1}+\mathrm{}+X_{[Nt]}^{\mu ^H,M(N)}}{N^H\sqrt{M(N)}}=B_H(t)$$

Proof: Using the generalization of Berry Esseen’s inequality to multidimensional variables , we get the convergence for finite-dimensional marginals.

To get the weak convergence, we cannot use directly  as in previous section, and we prove the tightness of the family of processes, by checking Billingsley’s criteria  (Theorem 15.6):

Denote by $`S_N(t):=\frac{X_{[Nt]}^{\mu ^H,1}+\mathrm{}+X_{[Nt]}^{\mu ^H,M(N)}}{\sqrt{M(N)}}`$.

Let $`1t_2tt_10`$, and $`k`$,

$`\begin{array}{cc}\hfill J_N(k,t_2,t,t_1)& :=E[|\frac{S_N(t_2)S_N(t)}{N^H}|^k|\frac{S_N(t)S_N(t_1)}{N^H}|^k]\hfill \\ & \frac{1}{N^{2kH}}E[S_N(t_2t)^{2k}]^{\frac{1}{2}}E[S_N(tt_1)^{2k}]^{\frac{1}{2}}\hfill \end{array}`$

$$k,E[S_N(t)^{2k}]=\frac{1}{M(N)^k}E[(X_{[Nt]}^{\mu ^H,1}+\mathrm{}+X_{[Nt]}^{\mu ^H,M(N)})^{2k}]$$

When we develope the polynomial with degree $`2k`$ inside the expectation, we notice that only the monomials of the type $`(X^1)^{2\alpha _1}\mathrm{}(X^M)^{2\alpha _M}`$ have a non null contribution, as far as the $`X^i`$’s are centered and independent. This set has cardinal $`O(M^k)`$.

Now, we find, by the same means than in Proposition 1, that

$$\alpha _11,E[X_N^{2\alpha _1}]=O(N^{2\alpha _1+(2H2)}).$$

This leads to $`E[S_N(t)^{2k}]=O((Nt)^{2k+(2H2)k})=O((Nt)^{2Hk}).`$

Hence, for some positive constant $`C`$,

$`J_N(k,t_2,t,t_1)C(t_2t)^{2Hk}(tt_1)^{2Hk}C(t_2t_1)^{2Hk}`$.

We choose $`k>\frac{1}{2H}`$ in order to satisfy Billingsley’s criteria.∎

We conjecture that this result remains true in the case of  which proposes an answer to the question asked at the end of .

We do the same for the case $`H=\frac{1}{2}`$:

###### Proposition 6.

For $`N`$ large enough,

$$E[(\frac{c_{\frac{1}{2}}|X_N^{\mu ^{\frac{1}{2}}}|}{\sqrt{N\mathrm{log}N}})^3]c_{\frac{1}{2}}\times \frac{\sqrt{2N}}{\mathrm{log}N}.$$

Proof: As in previous proposition, $`X_N=_{k=1}^N\epsilon _k`$, with:

_ $`\epsilon _n`$ are Bernoulli(1/2),

_ $`c_{\frac{1}{2}}^2Cov(\epsilon _k,\epsilon _l)=r(|kl|)`$, with $`r(n)_{\mathrm{}}\frac{1}{2n}`$.

Similarly as in previous proposition, $`E[({\displaystyle \frac{c_{\frac{1}{2}}|X_N|}{\sqrt{N\mathrm{log}N}}})^4]^{\frac{1}{2}}`$ is equivalent to:

$$\frac{c_{\frac{1}{2}}^2}{N\mathrm{log}N}(\frac{4!}{c_{\frac{1}{2}}^2}_0^N_0^{x_1}_0^{x_2}_0^{x_3}r((x_1x_2)+(x_3x_4))𝑑x_4𝑑x_3𝑑x_2𝑑x_1)^{\frac{1}{2}}$$

But,

$$_0^N_0^{x_1}_0^{x_2}_0^{x_3}\frac{1}{(x_1x_2)+(x_3x_4)}𝑑x_4𝑑x_3𝑑x_2𝑑x_1=\frac{N^3}{6}.$$

Now applying Berry-Esseen’s inequality, the error of the distribution function of the marginal at time 1, is dominated by $`1.3\times \frac{\sqrt{N}}{\mathrm{log}N\sqrt{M}}`$.

###### Corollary 2.

Let $`M`$ be a function on the integers such that $`\frac{M(N)}{N\mathrm{log}(N)^2}`$ tends to $`\mathrm{}`$,

$$^𝒟\underset{N\mathrm{}}{lim}c_{\frac{1}{2}}\frac{X_{[Nt]}^{\mu ^{\frac{1}{2}},1}+\mathrm{}+X_{[Nt]}^{\mu ^{\frac{1}{2}},M(N)}}{\sqrt{M(N)\times N\mathrm{log}(N)}}=B(t)$$

Finally, we treat the case $`H<\frac{1}{2}`$:

###### Proposition 7.

For $`H<1/2`$, for $`N`$ large enough,

$$E[(\frac{c_H|Y_N^{\mu ^H}|}{N^H})^3]D_HN^{\frac{1}{2}H},$$

with $`D_H=\sqrt{\frac{2H}{2H+1}}\times c_H.`$

Proof: The proof is more delicate than in Proposition 5, because the computation of the fourth moment makes appear compensations (as in the variance computation), and we have to treat them first ”$`p`$ by $`p`$”, before integrating against $`\mu ^H`$.

Assume $`i_4i_3i_2i_1`$, we get, as in Proposition 1, by successive conditionings:

$$E[\delta _{i_1}^p\delta _{i_2}^p\delta _{i_3}^p\delta _{i_4}^p]=r_p(i_1i_2)r_p(i_3i_4),$$

where $`r_p(0)=1`$ and $`n>0,r_p(n)=p(12p)^{n1}`$.

Now,

$$\begin{array}{cc}\hfill \underset{1i_kN}{}_{1k4}E[\delta _{i_1}^p\delta _{i_2}^p\delta _{i_3}^p\delta _{i_4}^p]& =12\underset{\mathrm{max}\{i_3,i_4\}\mathrm{min}\{i_1,i_2\}}{}E[\delta _{i_1}^p\delta _{i_2}^p\delta _{i_3}^p\delta _{i_4}^p]+O(N)\hfill \end{array}$$

$$=12\underset{\mathrm{max}\{i_3,i_4\}\mathrm{min}\{i_1,i_2\}}{}r_p(|i_1i_2|)r_p(|i_3i_4|)+O(N)$$

$$=12\underset{\mathrm{min}\{i_1,i_2\}=1}{\overset{N}{}}\left(\underset{\mathrm{max}\{i_3,i_4\}=1}{\overset{\mathrm{min}\{i_1,i_2\}}{}}r_p(|i_3i_4|)\right)r_p(|i_1i_2|)+O(N)$$

$$=12\underset{\mathrm{min}\{i_1,i_2\}=1}{\overset{N}{}}\left(\underset{\mathrm{max}\{i_3,i_4\}=1}{\overset{\mathrm{min}\{i_1,i_2\}}{}}\left(r_p(0)+2\underset{l=1}{\overset{\mathrm{max}\{i_3,i_4\}}{}}r_p(l)\right)\right)r_p(|i_1i_2|)+O(N)$$

$$=12\underset{\mathrm{min}\{i_1,i_2\}=1}{\overset{N}{}}\left(\underset{\mathrm{max}\{i_3,i_4\}=1}{\overset{\mathrm{min}\{i_1,i_2\}}{}}(12p)^{\mathrm{max}\{i_3,i_4\}}\right)r_p(|i_1i_2|)+O(N)$$

$$=12\underset{\mathrm{min}\{i_1,i_2\}=1}{\overset{N}{}}\left(\underset{k=1}{\overset{\mathrm{min}\{i_1,i_2\}}{}}(12p)^k\right)\left(r_p(0)+2\underset{l=1}{\overset{N\mathrm{min}\{i_1,i_2\}}{}}r_p(l)\right)+O(N)$$

$$=12\underset{\mathrm{min}\{i_1,i_2\}=1}{\overset{N}{}}\left(\underset{k=1}{\overset{\mathrm{min}\{i_1,i_2\}}{}}(12p)^k\right)(12p)^{N\mathrm{min}\{i_1,i_2\}}+O(N)$$

$$=12\underset{1klN}{}(12p)^{N(lk)}+O(N)$$

(The contribution of the exceptional situations $`\mathrm{max}\{i_3,i_4\}=\mathrm{min}\{i_1,i_2\}`$ is estimated by $`O(N)`$, because of the same compensations as the ones described in the above equalities.)

We now use,

$$\begin{array}{cc}\hfill _0^1(12p)^n𝑑\mu ^H(p)& =2\underset{k>n}{}r(k)\hfill \\ & \frac{1}{c_H^2}\frac{2H}{n^{12H}}\hfill \end{array}$$

We deduce,

$$\frac{c_H^2}{N^{2H}}(\underset{1i_kN}{}_{1k4}E[\delta _{i_1}\delta _{i_2}\delta _{i_3}\delta _{i_4}])^{\frac{1}{2}}$$

$$\frac{c_H^2}{N^{2H}}\left(\frac{2H}{c_H^2}_0^N_0^{x_2}\frac{dx_1}{(N(x_2x_1))^{12H}}\right)^{\frac{1}{2}}$$

$$=c_H\times \sqrt{\frac{2H}{2H+1}}N^{\frac{1}{2}H}.$$

Applying Berry-Esseen’s inequality, the error of the distribution function of the marginal at time 1, is dominated by $`0.65\times D_H\frac{N^{\frac{1}{2}H}}{\sqrt{M}}`$.

###### Corollary 3.

Let $`M`$ be a function on the integers such that $`\frac{M(N)}{N^{12H}}`$ tends to $`\mathrm{}`$,

$$^𝒟\underset{N\mathrm{}}{lim}c_H\frac{X_{[Nt]}^{\mu ^H,1}+\mathrm{}+X_{[Nt]}^{\mu ^H,M(N)}}{N^H\sqrt{M(N)}}=B_H(t)$$

As a general conclusion, we can say that, using $`\mu _H`$, the number of computations we have to make, in order to get a $`N`$ steps trajectory, is a constant times $`N^{32H}`$ for $`\frac{1}{2}<H`$, and $`N^{22H}`$ for $`0<H<\frac{1}{2}`$. In any case, it is a power of $`N`$ between 1 and 2.

The constant in factor, is of big importance for practical simulation. We took the best constant in the Berry-Esseen’s inequality we could find in the literature, even if it is bigger than the best constant possible conjectured by Esseen, i.e. $`\frac{3+\sqrt{10}}{6\sqrt{2\pi }}0.41`$ that would gain in time a squared factor equal to 2.5 (see  for a nice discussion on this subject).

But more important is to note that the constant can be considerably ameliorated by using other measures than $`\mu _H`$, providing smaller $`c_H`$’s. It is the case for the real-indexed sequence of probabilities $`(\mu _{H,k})_{k>0}`$, defined below:

_ $`H[\frac{1}{2},1],`$

$$d\mu _{H,k}(p):=2^{k+12H}\frac{\mathrm{\Gamma }(k+22H)}{\mathrm{\Gamma }(k)\mathrm{\Gamma }(22H)}(p\frac{1}{2})^{k1}(1p)^{12H}1_{[\frac{1}{2},1]}(p)dp$$

which is just the law of $`\frac{1+B(k,22H)}{2}`$ and coincides with $`\mu _H`$ for $`k=1`$ (where $`B(a,b)`$ denotes the Beta variable with parameters $`a`$ and $`b`$).

_ $`H]0,\frac{1}{2}[,`$

$$d\mu _{H,k}(p):=2^{k2H}\frac{\mathrm{\Gamma }(k+12H)}{\mathrm{\Gamma }(k)\mathrm{\Gamma }(12H)}(\frac{1}{2}p)^{k1}p^{2H}1_{[0,\frac{1}{2}]}(p)dp$$

which is just the law of $`\frac{B(12H,k)}{2}`$ and coincides with $`\mu _H`$ for $`k=1`$.

_ For $`H>\frac{1}{2}`$, we obtain $`c_{H,k}=\sqrt{\frac{H(2H1)\mathrm{\Gamma }(k)}{\mathrm{\Gamma }(k+22H)}}\frac{\sqrt{H(2H1)}}{k^{1H}}`$.

So that, using $`\mu _{H,k}`$ for large $`k`$, yields an error equivalent to:

$$0.65\sqrt{\frac{6H(2H1)^2}{(H+1)(2H+1)}}\times (\frac{N}{k})^{1H}\frac{1}{\sqrt{M}}$$

_ For $`H=\frac{1}{2}`$, we obtain $`c_{\frac{1}{2},k}=\frac{1}{\sqrt{2k}}`$.

The error is equivalent to:

$$\frac{0.65}{\mathrm{log}N}\times (\frac{N}{k})^{\frac{1}{2}}\frac{1}{\sqrt{M}}$$

_ For $`H<\frac{1}{2}`$, we obtain $`c_{H,k}=\sqrt{\frac{2H\mathrm{\Gamma }(k)}{\mathrm{\Gamma }(k+12H)}}\frac{\sqrt{2H}}{k^{\frac{1}{2}H}}`$.

The error is equivalent to:

$$0.65\sqrt{\frac{4H^2}{2H+1}}\times (\frac{N}{k})^{\frac{1}{2}H}\frac{1}{\sqrt{M}}$$

Using $`\mu _{H,k}`$ instead of $`\mu _H`$, allows to gain for $`M`$ a factor $`k^{22H}`$ (resp. $`k^{12H}`$) for $`H>\frac{1}{2}`$ (resp. for $`H<\frac{1}{2}`$). Loosely speaking, it erases the noise generated by the parameters between $`\frac{1}{2}`$ and a fixed constant smaller than 1. The trouble making $`k`$ increase, is that it damages the value of the covariance of $`X_N`$ (resp. $`Y_N`$), but we can allow any $`k(N)=o(N)`$. So that we obtain an algorithm with any $`M(N)=1/o(1)`$ number of trajectories. As a result, our algorithm requires a number of computations of the order $`N/o(1)`$, for any $`o(1)`$. Moreover, we have only $`M(N)`$ real datas to keep in memory along the whole procedure, corresponding to the coin-tossed parameters of the walks, and $`M(N)`$ integers (-1 or 1), giving the last moves of the walks.

We want to give now a second family of measures $`(\mu _{H,k}^{})_{k>0}`$:

_ For $`H>\frac{1}{2}`$, $`\mu _{H,k}^{}`$ is the law of $`1\frac{(1U^{\frac{1}{k}})^{\frac{1}{22H}}}{2}`$.

An easy computation gives $`c_{H,k}^{}=\frac{c_H}{\sqrt{k}}`$ (error: $`0.65\times D_H\frac{N^{1H}}{\sqrt{kM}}`$).

_ For $`H<\frac{1}{2}`$, $`\mu _{H,k}^{}`$ is the law of $`\frac{(1U^{\frac{1}{k}})^{\frac{1}{12H}}}{2}`$.

Again, an easy computation gives $`c_{H,k}^{}=\frac{c_H}{\sqrt{k}}`$ (error: $`0.65\times D_H\frac{N^{\frac{1}{2}H}}{\sqrt{kM}}`$).

The advantage of this family is obviously the easy simulation it provides. The error is estimated by a term containing $`\frac{1}{\sqrt{k}}`$, that seems to be better than the last one, but the damages on the variances grow faster than for $`\mu _{H,k}`$. Actually the scale $`\sqrt{k}`$ corresponds, in the previous family, to the scale $`k^{1H}`$ (resp. $`k^{\frac{1}{2}H}`$) for $`H>\frac{1}{2}`$ (resp. $`H<\frac{1}{2}`$). The drawback of this family, is that the theoretical computations of the variance are not very explicit.

We illustrate our results by three graphs corresponding to three different parameters of $`H`$, with $`N=1000`$, and a theoretical error smaller than $`10\%`$ (we indicate the time it takes for Matlab to draw a graph):

_ For $`H=0.25`$, we take $`M=200`$ and $`k=1`$ (15 seconds) (Fig. 1).

_ For $`H=0.5`$, we can see that the convergence in Theorem 2 is the slowest one, and we will use the simple random walk to simulate it. (Fig. 2)

_ For $`H=0.75`$, we take $`M=400`$ and $`k=0.5`$ (25 seconds) (Fig. 3).

Remark: The interest of taking large $`k`$ appears when $`N`$ is very large, especially if we want to preserve a small error on the whole trajectory. In the case $`H>\frac{1}{2}`$, we were even obliged to take $`k`$ smaller than one.

We remark the different behaviours of the trajectories: let us recall that the Hausdorff dimension of the trajectories are a.s. equal to $`2H`$, and we notice that the variances of the process between 0 and 1 become larger when $`H`$ decreases.

As we noticed just above, we may be limited by the fidelity of the covariance of our process. In this spirit, it is quite interesting to note that the autocovariance function of the Gaussian noise $`G_H(j)`$ of the introduction is, up to a shift, the sequence of moments of a probability measure on $`[0,1]`$. I first remarked it by checking the conditions of Hausdorff theorem ( p.226), but Marc Yor gave me kindly, the method to get the explicit density of this measure, and I present it here. This brings a third (the last !) family of probability measures:

###### Proposition 8.

Let $`H]\frac{1}{2},1[`$. Consider the family of probability measures $`(\nu _{H,k})_{k>0}`$ on $`[\frac{1}{2},1]`$, with density $`C(H,k)\times (1p)^2(2p1)^{k1}(\mathrm{ln}(\frac{1}{2p1}))^{12H}`$, and

$`C(H,k):=\frac{16H(2H1)}{\mathrm{\Gamma }(22H)}\times ((k+2)^{2H}2(k+1)^{2H}+k^{2H})^1`$

$$n0,_0^1(2p1)^n𝑑\nu _{H,k}(p)=\frac{(n+k+2)^{2H}2(n+k+1)^{2H}+(n+k)^{2H}}{(k+2)^{2H}2(k+1)^{2H}+k^{2H}}$$

Proof: $`n1`$,

$$\frac{1}{2}((n+1)^{2H}2n^{2H}+(n1)^{2H})=H_0^1(n+t)^{2H1}(n+t1)^{2H1}dt$$

$$=H(2H1)_0^1_0^1(n+t+s1)^{2H2}𝑑s𝑑t$$

$$=\frac{H(2H1)}{\mathrm{\Gamma }(22H)}_0^1_0^1_0^+\mathrm{}e^{(n+t+s1)u}u^{12H}𝑑u𝑑s𝑑t$$

$$=\frac{H(2H1)}{\mathrm{\Gamma }(22H)}_0^+\mathrm{}e^{nu}(1e^u)^2e^uu^{12H}𝑑u$$

$$=\frac{H(2H1)}{\mathrm{\Gamma }(22H)}_0^1x^n(\frac{1x}{x})^2(\mathrm{ln}\frac{1}{x})^{12H}𝑑x$$

The end of the proof is straightforward, by change of variable.∎

Note: The normalization corresponding to $`\nu _{H,k}`$ is:

$`c_{H,k}^{\prime \prime }=(\frac{(k+2)^{2H}2(k+1)^{2H}+k^{2H}}{2})^{\frac{1}{2}}\frac{\sqrt{H(2H1)}}{k^{1H}}c_{H,k}`$.

The same can be done for $`H]0,\frac{1}{2}[`$:

###### Proposition 9.

Let $`H]0,\frac{1}{2}[`$. Consider the family of probability measures $`(\nu _{H,k})_{k>0}`$ on $`[0,\frac{1}{2}]`$, with density $`C(H,k)\times p(2p1)^{k1}(\mathrm{ln}(\frac{1}{12p}))^{12H}`$, and

$`C(H,k):=\frac{8H}{\mathrm{\Gamma }(12H)}\times ((k+1)^{2H}k^{2H})^1`$

$$n1,_0^1p(2p1)^{n1}𝑑\nu _{H,k}(p)=\frac{(n+k+1)^{2H}2(n+k)^{2H}+(n+k1)^{2H}}{2((k+1)^{2H}k^{2H})}$$

Proof: $`n1`$,

$$\frac{1}{2}((n+1)^{2H}2n^{2H}+(n1)^{2H})=H_0^1(n+t)^{2H1}(n+t1)^{2H1}dt$$

$$=\frac{H}{\mathrm{\Gamma }(12H)}_0^1_0^+\mathrm{}(e^{(n+t)u}e^{(n+t1)u})u^{2H}𝑑u𝑑t$$

$$=\frac{H}{\mathrm{\Gamma }(12H)}_0^+\mathrm{}e^{(n+1)u}(e^u1)^2u^{12H}𝑑u$$

$$=\frac{H}{\mathrm{\Gamma }(12H)}_0^1x^n(\frac{1x}{x})^2(\mathrm{ln}\frac{1}{x})^{12H}𝑑x$$

$$=\frac{8H}{\mathrm{\Gamma }(12H)}_0^{\frac{1}{2}}(12p)^{n2}p^2(\mathrm{ln}\frac{1}{12p})^{12H}𝑑p$$

We find the normalization, using $`_{n1}(n+k+1)^{2H}2(n+k)^{2H}+(n+k1)^{2H}=k^{2H}(k+1)^{2H}`$, and the relation $`1=2_{n1}p(12p)^{n1}`$.∎

Note: The normalization corresponding to $`\nu _{H,k}`$ is:

$`c_{H,k}^{\prime \prime }=((k+1)^{2H}k^{2H})^{\frac{1}{2}}\frac{\sqrt{2H}}{k^{\frac{1}{2}H}}c_{H,k}`$.

These measures $`\nu _{H,k}`$ have an interest, as far as $`X_N^{\nu _H}`$ (resp. $`Y_N^{\nu _H}`$) have explicit covariances matching quite well with the covariances of the fractional Brownian motion. On the other hand, they don’t seem very easy to simulate.

Fig. 1: H=0.25

Fig. 2: H=0.5

Fig. 3: H=0.75

Aknowledgements: First of all, I would like to thank Erick Herbin for asking me the question about the relation between fractional Brownian motion and correlated random walks. It is a pleasure for me to thank also Marc Yor who helped me finding $`\nu _H`$, and Zhan Shi for fruitful conversations.