Two-species diffusion-annihilation process on the fully-connected lattice: probability distributions and extreme value statistics
Loïc Turban














|#⟩1|#1⟩
⟨#|1⟨#1|
⟨#|1⟩#2⟨#1|#2⟩
0pt minus 1fil


#1#1

1

     
     
  Keywords: #1


    Δ𝐄𝐀𝐁𝐈𝒰𝒱𝒥𝒦ℒℛ𝒯𝒮𝒩𝒟𝔗𝔔𝔙𝔖𝔓Δs.l.=s.l.≃D_seg1Two-species diffusion-annihilation process]Two-species diffusion-annihilation process 
        on the fully-connected lattice: probability distributions and extreme value statisticsLaboratoire de Physique et Chimie Théoriques, Université de Lorraine–CNRS (UMR7019),
        Vandœuvre lès Nancy Cedex, F-54506, Franceloic.turban@univ-lorraine.fr We study the two-species diffusion-annihilation process, A+B→ Ø,
        on the fully-connected lattice. Probability distributions for the number of particles and 
        the reaction time are obtained for a finite-size system using a master equation approach. 
        Mean values and variances are deduced from generating functions. When the reaction is far 
        from complete, i.e. for a large number of particles of each species, mean-field theory is exact
        and the fluctuations are Gaussian. 
        In the scaling limit the reaction time displays extreme-value statistics in the vicinity of 
        the absorbing states. A generalized Gumbel distribution is obtained for unequal initial densities, 
        ρ_A>ρ_B. For equal or almost equal initial densities, 
        ρ_A≃ρ_B, the fluctuations of the reaction time near the absorbing state 
        are governed by a probability density involving derivatives of ϑ_4, 
        the Jacobi theta function.
        J. Phys. A: Mathematical and Theoretical
        
        § INTRODUCTION
        
        
        In the field of non-equilibrium statistical mechanics, reaction-diffusion processes 
        offer the possibility to study the effects of fluctuations on conceptually very 
        simple model systems like the single-species or the two-species annihilation 
        processes <cit.>.
        
        In a standard mean-field approximation <cit.>, the bimolecular reaction 
        A+B→ Ø displays a t^-1 asymptotic decay of the particle densities 
        for equal initial values, ρ_A(0)=ρ_B(0). For unequal densities, ρ_A>ρ_B, 
        the approach to the 
        absorbing state, ρ_B=0 and ρ_A=ρ_A(0)-ρ_B(0), is exponential.
        The mean-field approximation assumes that the system remains homogeneous and ignores 
        the effect of spatial correlations in the distribution of reactants, thus giving a lower bound to
        the actual particle densities <cit.>.
        
        The relevance in low dimensions of initial concentration fluctuations was pointed out
        by Ovchinnikov and Zeldovich <cit.> who found a t^-3/4 decay in dimension D=3
        for equal initial densities.
        This result was soon generalized and a t^-D/4 
        decay was proposed for ρ_A(0)=ρ_B(0) on the basis of numerical simulations, approximate 
        analytical approaches and scaling 
        arguments <cit.> [
        Note that for initially separated reactants the kinetics is inhomgeneous and governed 
        by the reaction in a growing domain around the interface <cit.>.]. 
        The validity of this asymptotic behaviour was later confirmed by establishing
        rigorous bounds on the particle density <cit.> and 
        through a renormalization group study <cit.>. 
        
        The slowing down of the process is due to the segregation of A and B particles
        into A-rich and B-rich domains, at the scale of the diffusion 
        length <cit.>.
        The segregation is a consequence of the initial fluctuations of the densities around their mean values.
        At long time the reaction is efficient only at the interface between the domains and thus slows down.
        This effect is relevant below the segregation dimension =4 
        at which the t^-1 homogeneous mean-field decay is recovered. 
        When generalized to q species <cit.> the problem has a segregation dimension
        (q)=4/(q-1)≥2 <cit.>.
        
        When ρ_A(0)>ρ_B(0) the density of the minority reactant behaves asymptotically as 
        ρ_B(t)∼^-λ_DG_D(t) ,
        
        with <cit.>
        G_D(t)={[   √(t) ,      D=1; t/ln t ,      D=2;      t ,      D=3 ].
        
        Note that the upper critical dimension, 
        as for the single-species process, is D_c=2 <cit.>. The 
        t^-D/4 behaviour can be actually obtained using mean-field rate equations, provided the 
        inhomogeneity of the system is taken into account. Although there is no qualitative change
        at D_c for equal initial densities, the upper critical dimension signals itself
        via logarithmic corrections at D_c and a stretched exponential decay below 
        D_c for unequal initial densities.
        
        The two-species annihilation process has potential applications in different domains. It can be 
        used to model particle-antiparticle 
        annihilation in the early universe <cit.>, the kinetics of
        bimolecular chemical reactions <cit.> or 
        electron-hole recombination in irradiated semiconductors <cit.>.
        
        The aim of the present work is to study analytically the kinetics 
        of the two-species reaction-diffusion process 
        on the fully-connected lattice with an emphasis on probability distributions. 
        This is a continuation of previous work on the single-species 
        process <cit.>. Since the lattice with N sites can only be embedded in a (N-1)-dimensional 
        space, taking the thermodynamic limit requires an infinite-dimensional space and one expects mean-field 
        behaviour. Our purpose is to obtain exact results for the particle density and the 
        reaction time in finite-size systems and to study the 
        extreme-value statistics of the reaction time, in the vicinity of the absorbing state, 
        for both equal or unequal initial densities of the reactants.
        
        The paper is organized as follows. In section 2 we present the model, its mean-field solution when 
        homogeneity is assumed and give a brief description of our results. In section 3
        we study the statistics of the number of particles surviving at a given time, first on a finite system 
        and then in the scaling limit. Section 4 is devoted to a similar 
        study of the reaction time, i.e. the time 
        needed to have a given number of particles remaining. This is followed by the conclusion in section 5. 
        Details of the calculations are given in six appendices.
        
        
        
        § MODEL, MEAN FIELD AND MAIN RESULTS
        
        
         §.§ Model
        
        
        
        
        We consider the two-species reaction-diffusion process, A+B→ Ø, 
        on a fully connected lattice with N sites. 
        Let s_A and s_B be the number of particles of each type with s_A≥ s_B 
        and at most one particle per site. In the following we shall use the variables
        
        s=s_A+s_B/2 ,     d=s_A-s_B/2 .
        
        Thus d=0 when the initial densities are equal and s=d when the reaction is complete. 
        A dictionary giving the relations with standard notations is given in table <ref>.
        
        The system evolves in time through random sequential updates. 
        An update consists of one or two steps. A first site i is selected at random among the N.  
        When this site is occupied by a particle of type A (B) a second site j is randomly 
        selected among the N. If the destination site is occupied by a particle of type B (A), 
        the two particles annihilate and s→ s-1 . In all other cases s is unchanged. 
        d is always conserved. At each update the time t is incremented 
        by 1/N so that t=k/N where k is the number of updates. Note that first selecting a site
        instead of a particle is vital to keep a constant time increment.
        
        The probabilities for the different events are the following:
        
          * s→ s'=s-1, with probability:
        
        2 s^2-d^2/N^2=
        s_A/N_i=A×s_B/N_j=B 
        + s_B/N_i=B×s_A/N_j=A
          * s→ s'=s, with probability:
        
        1-2 s^2-d^2/N^2=1-2s/N_i=∅ + s_A/N_i=A×1-s_B/N_j≠B + s_B/N_i=B×1-s_A/N_j≠A
        Note that, contrary to what occurs on finite-dimensional lattices, 
        the initial distribution of the two species does not matter for the fully-connected lattice.
        
        
        
         §.§ Mean field solution
        
        In the following we always assume that in the initial state, at t=0, all the sites are occupied, 
        s_A+s_B=N and we neglect the effect of spatial fluctuations (D>). The initial value of s is then n=N/2 according to sd. 
        In the scaling limit (s.l.), when N and n→∞, we introduce the scaled variables
         
        xs/n=ρ_A+ρ_B ,     
        yd/n=ρ_A-ρ_B .
        
        where ρ_A and ρ_B are the particle densities.
        Note that x is the fraction of occupied sites at t and y 
        is a constant giving the asymptotic value of this fraction when t→∞.
        
        After a small number of updates, Δ k, according to probrea 
        the mean value of s is changed by 
        Δ s=-2 s^2-d^2/N^2Δ k 
        where on the right the fluctuations of s 
        around its mean value are neglected. In the scaling limit,
        Δ s→ n dx, Δ k→ N dt=2n dt, yielding
        dx/dt=-s^2-d^2/n^2=-(x^2-y^2) ,
        
        so that
        dx/x^2-y^2=1/2y(dx/x-y-dx/x+y)=-dt .
        
        Finally the solution satisfying the initial condition, x=1 when t=0, is given by:
        
        x=y 1+y+(1-y) ^-2yt/1+y-(1-y) ^-2yt .
        
        This yields
        ρ_A   =s_A/N=x+y/2=y(1+y)/1+y-(1-y) ^-2yt ,
        ρ_B   =s_B/N=x-y/2=y(1-y) ^-2yt/1+y-(1-y) ^-2yt ,
        
        for the densities of the two species. When y>0 the approach to the asymptotic values, 
        ρ_A=y and ρ_B=0, is exponential. When y→0 an algebraic decay is obtained:
        ρ_A=ρ_B=1/2(t+1) .
        
        
        
        
        The time needed to reach a given value of x is
        
        t=1/2yln[(1-y)(x+y)/(1+y)(x-y)] ,
        
        leading to
        
        t=1/x-1 ,
        
        when y→0.
        
        
        
         §.§ Main results
        
        
        
        The evolution of the system is illustrated in figures <ref> and <ref> 
        giving respectively the probability distribution S_N(s,k) of s at different reaction times and 
        the probability distribution T_N(s,k) of the reaction time at different values of s. 
        
        
        As expected for an infinite-dimensional system the mean value of s at time t 
        is in agreement with mean-field theory. Asymptotically, it decays as 
        t^-1 when s_A=s_B and approaches its asymptotic value exponentially when s_A>s_B. 
        The mean value and the variance both scale as n and the fluctuations of s are Gaussian.
        
        For the statistics of the reaction time t three different 
        regimes are observed with the following results in the scaling limit: 
        
        
          * x>y≥0: The reaction is far from complete with s_A=(n), s_B=(n) and s_A≥ s_B. 
        The mean values of the reaction time are the mean-field ones:
        t_N=1/x-1 ,   y=0 ;    t_N=1/2yln[(1-y)(x+y)/(1+y)(x-y)] ,   y>0 .
        
        The fluctuations are weak, the variance scaling as n^-1:
        Δ t_N^2   =χ_y(x)/n ;    χ_0(x)=1/3x^3-1/2x+1/6 ,   y=0 ;
        χ_y(x)   =1/2y^2(x/x^2-y^2-1/1-y^2)
        -1/4y(1+1/y^2)
        ln[(1-y)(x+y)/(1+y)(x-y)] ,   y>0 .
        
        The probability density is Gaussian:
        (x,θ)=^-θ^2/[2χ_y(x)]/√(2πχ_y(x)) ,  θ=n^1/2(t-t_N) .
        
          * x=y>0, s=ny+u:  
        The reaction is close to completion with unequal numbers of  
        particles (s_A-s_B=(n) and u=s_B=(1)). The reaction time scales logarithmically with nt_N=1/2y[ln(2ny1-y/1+y)
        +γ-H_u] ,
        
        where H_u=∑_j=1^u1/j is a harmonic number and  γ=0.577215665… is the Euler constant. 
        The variance is independent of n:
        Δ t_N^2=1/4y^2[ζ(2)-H_u^(2)] .
        
        Here H_u^(2) is a generalized harmonic number such that H_l^(m)=∑_j=1^l1/j^m
        and ζ(2)=π^2/6.
        The system displays extreme value statistics. The fluctuations are governed by a generalized Gumbel 
        distribution <cit.>, indexed by u:
        '(u,θ')=1/u!exp[-(u+1)(θ'+γ-H_u)
        -^-(θ'+γ-H_u)] ,  θ'=2y(t-t_N) .
        
          * x=y=0, s≥ d: 
         The reaction is close to completion with s_A=(1)≥ s_B≥0.
        The reaction time grows as nt_N=n(ζ(2)-H_s^(2)) ,   d=0 ;    t_N=n/2d(H_s+d-H_s-d) ,   d>0 ,
        
        and the variance as n^2Δ t_N^2   =n^2(ζ(4)-H_s^(4)) ,   d=0 ;
        Δ t_N^2   =n^2/4d^2[2ζ(2)-H_s-d^(2)-H_s+d^(2)+1/d(H_s-d-H_s+d)]
         ,   d>0 ,
        
        where ζ(4)=π^4/90. The fluctuations of the reaction time are even stronger 
        and  governed now by derivatives of the Jacobi theta function ϑ_4:
        ”(s,θ”)=^d^2θ”/(s-d)!(s+d)!∏_m=0^s(m^2+d/dθ”) ϑ_4(0,^-θ”) ,  θ”=t/n .
        
        
        § NUMBER OF SURVIVING PARTICLES AT A GIVEN TIME
        
        In this section we study the probability distribution S_N(s,k) giving the probability to have s+d 
        particles of type A and s-d particles of type B remaining after k updates. As above we assume 
        that the N sites are initially occupied, s_A+s_B=N.
        
        
        
         §.§ Master equation
        
        According to probrea and probnorea the master equation governing the evolution of 
        the system takes the following form
        
        S_N(s,k)=(1-2 s^2-d^2/N^2)S_N(s,k-1)
        +2 (s+1)^2-d^2/N^2 S_N(s+1,k-1) ,   s=d,…,n ,
        
        with the boundary condition S_N(s>n,k)=0 and the initial condition S_N(s,0)=δ_s,n.
        In master the first (second) term on the right gives the probability to 
        be in a state with 2s particles (2s+2 particles) after k-1 updates and to remain in 
        this state (to have two particles annihilating) at the kth update.
        
        
        
         §.§ Eigenvalue problem
        
        Let us define the column state vector |S_N(k)⟩ with components S_N(s,k),   s=d,…,n,
        the master equation master can be written in matrix form as 
        |S_N(k)⟩=𝖳|S_N(k-1)⟩ where the transition matrix 𝖳 is given by:
        𝖳=([                   1   2 (d+1)^2-d^2/N^2                   0                   0                   0                   0;                   0 1-2 (d+1)^2-d^2/N^2   2 (d+2)^2-d^2/N^2                   0                   0                   0;                                                           ⋱                   ⋱                                        ;                   0                   0                   0     1-2 s^2-d^2/N^2   2 (s+1)^2-d^2/N^2                   0;                                                                                                   ⋱                   ⋱;                   0                   0                   0                   0                   0     1-2 n^2-d^2/N^2;                     ]) .
        
        The eigenvalue equation 𝖳|v^(r)⟩=λ_r|v^(r)⟩ leads to the linear system
        (1-2 s^2-d^2/N^2-λ_r)v_s^(r)+2 (s+1)^2-d^2/N^2 v_s+1^(r)=0 ,
           s=d,…,n ,
        
        with v_n+1^(r)=0. It is easy to verify that 
        λ_r=1-2 r^2-d^2/N^2 ,     
        v_s^(r)={[ (-1)^r-s
         v_r^(r)∏_j=1^r-s(s+j)^2-d^2/r^2-(s+j-1)^2                                                when                                                 s<r;                                                   0                                                when                                                 s>r ].,
        
        solves the eigenvalue problem lin. The solution involves the repeated use of the recursion relation
        
        v_s^(r)=-(s+1)^2-d^2/r^2-s^2 v_s+1^(r) ,
        
        which follows from lin.
        The value of v_r^(r), which remains free, will be used 
        to satisfy the initial condition.
        
        
        
         §.§ Probability distribution S_N(s,k)
        
        We look for the initial state vector under the form |S_N(0)⟩=∑_r=d^n|v^(r)⟩ 
        which leads to the condition
        
        S_N(s,0)=∑_r=s^nv_s^(r)=δ_s,n ,     s=d,…,n ,
        
        for the components. From the values of v_r^(r) with r=n,…,n-3 (see appendix A) 
        we can infer that the general expression reads:
        
        v_r^(r)=(n-d)!(n+d)!/(n-r)!(n+r)!2r r-d ,     r=d,…,n .
        
        Then, according to sol, one obtains:
        
        v_s^(r)=(-1)^r-s2r(n-d)!(n+d)!(r+s-1)!/(n-r)!(n+r)!(s-d)!(s+d)!(r-s)! ,   s≤ r .
        
        After k updates the state vector |S_N(k)⟩ is given by
        𝖳^k|S_N(0)⟩=∑_r=d^n𝖳^k|v^(r)⟩
        =∑_r=d^nλ_r^k|v^(r)⟩ ,
        
        which, according to sol and vsr, gives
         S_N(s,k)=∑_r=s^nλ_r^kv_s^(r)
        =∑_r=s^n
        (-1)^r-s2r(n-d)!(n+d)!(r+s-1)!/(n-r)!(n+r)!(s-d)!(s+d)!(r-s)!(1-2 r^2-d^2/N^2)^k
        
        for the components.
        
        
        
         §.§ Mean value and variance when d=0
        
        
        
        Let us define the generating function
        _N(w,k)=∑_s=d^nsw^sS_N(s,k)=∑_r=d^nr(n-d)!(n+d)!/(n-r)!(n+r)!(1-2 r^2-d^2/N^2)^kΩ_r,d(w)
        
        where
        Ω_r,d(w)=∑_s=d^r(-1)^r-sr+s s-dr+d s+d2sw^s/r+s .
        
        In appendix B we show that when d=0Ω_r,0(1)=2
        
        and
        .dΩ_r,0/d w|_w=1=2r^2 ,
        
        which allows us to evaluate the mean value of ss_N(k)=_N(1,k)=∑_r=1^n2r∏_j=1^rn-j+1/n+j(1-2 r^2/N^2)^k ,   d=0 ,
        
        and its mean-square value:
        s_N^2(k)=.∂_N/∂ w|_w=1=∑_r=1^n2r^3∏_j=1^rn-j+1/n+j(1-2 r^2/N^2)^k ,   d=0 .
        
        In the scaling limit (N,n,k→∞, t=k/N) studied in appendix C, one obtains:
        s_N(t)/n1/t+1-1/6n(1-1/(t+1)^3) ,
          s_N^2(t)/n^21/(t+1)^2-t/2n(t+1)^4 ,   d=0 .
        
        In these expressions we kept the sub-leading contributions since the leading ones 
        vanish in the variance given by:
        Δ s_N^2(t)/n1/3(t+1)-1/2(t+1)^3+1/6(t+1)^4 ,
             d=0 .
        
        Thus the fluctuations are small and s is self-averaging. A comparison with finite-size data
        is shown in figure <ref>.
        We were not able to evaluate Ω_r,d(w) when d>0. This case is treated
        directly in the scaling limit in the next section.
        
        
        
         §.§ Scaling limit when d≥0
        
        
        
        Let us assume that in the scaling limit, for any value of d, s_N 
        and Δ s_N^2 are both growing as n as 
        in snt and dsn2t for d=0. This suggests the introduction, 
        besides the time variable t=k/(2n) and the density 
        y=d/n, of the scaled and centered variable
        σ(s,k)=s-s_N(k)/n^1/2 .
        
        Furthermore let us write the unknown mean value x as
        x=s_N/n=g_y(t) ,
        
        and define the probability density:
        (σ,t) n^1/2S_N(s,k) .
        
        Starting  from the master equation master with S_N replaced by , 
        a Taylor expansion of the right-hand-side
        up to second order in s and k, when re-expressed in terms of the scaled variables, 
        takes the form of an expansion in powers of n^-1/2. The terms independent of n 
        cancel. The terms of order n^-1/2 leads to the differential equation
        dg_y/dt=-(g_y^2-y^2) ,
        
        which is the mean-field equation edif-1 so that, according to xmf,
        x=g_y(t)=y 1+y+(1-y) ^-2yt/1+y-(1-y) ^-2yt ,
        
        in agreement with snt when y→0.
        To the next order, n^-1, one obtains the following partial differential equation:
        ∂/∂ t=-1/2[dg_y/dt
        +1/2(dg_y/dt)^2]∂^2/∂σ^2
        +2g_y(+σ∂/∂σ) .
        
        
        Introducing the reduced variance κ_y(t)=Δ s_N^2/n and assuming that 
         depends on t only through κ_y, the partial differential equation dsdt 
        can be rewritten in the following form:
        (∂/∂κ_y-1/2∂^2/∂σ^2)
        dκ_y/dt   =-1/2[dκ_y/dt+4g_yκ_y+dg_y/dt+1/2(dg_y/dt)^2]∂^2/∂σ^2
                                                  
        +2g_y(+σ∂/∂σ+κ_y∂^2/∂σ^2).
        
        The left-hand-side and the last bracket on the right-hand-side are both vanishing 
        for the Gaussian density with variance κ_y(t), thus we have
        (σ,t)=^-σ^2/[2κ_y(t)]/√(2πκ_y(t)) ,
        
        when, according to dsdkap, κ_y(t) satisfies the first-order differential equation:
        dκ_y/dt+4g_yκ_y+dg_y/dt+1/2(dg_y/dt)^2=0 .
        
        The solution is discussed in appendix D and reads:
        κ_y(t)Δ s_N^2/n=[(1+y)^2
        ^2yt-(1-y)^2^-2yt-4y(1+t-y^4t)]
        2y(1-y^2)^-4yt/[1+y-(1-y)^-2yt]^4 .
        
        In the limit y→0dsn2t is recovered. The mean value in xmean and the variance are 
        compared to finite-size data in figure <ref>.
        
        
        
        
        The scaling behaviour of the probability distribution S_N(s,k) at different times 
        is shown in figure <ref> for ρ_A=ρ_B and figure <ref> for 
        ρ_A-ρ_B=1/8.
        
        
        
        
        
        
        § TIME REQUIRED TO REACH A GIVEN NUMBER OF SURVIVING PARTICLES
        
        
         §.§ Probability distribution
        
        
        This section is dedicated to the study of T_N(s,k), the probability distribution for the 
        number of updates k needed to reach for the first time a total number 
        of surviving particles s_A+s_B=2s, 
        as shown in figure <ref>(a). This probability is related to S_N(s+1,k-1)
        through:
        
        T_N(s,k)=S_N(s+1,k-1)× 2(s+1)^2-d^2/N^2 .
        
        It is given by the product of the probability to be in a state with s'=s+1 
        after k-1 updates by the probability of the transition 
        s'=s+1→ s”=s at the next update.
        Making use of snsk one obtains:
        
        T_N(s,k)=(n-d)!(n+d)!/n^2(s-d)!(s+d)!∑_r=s+1^n
        (-1)^r-s-1r(r+s)!/(n-r)!(n+r)!(r-s-1)!(1-2 r^2-d^2/N^2)^k-1 .
        
        The evolution with k is governed by the master equation
        
        T_N(s,k)=[1-2 (s+1)^2-d^2/N^2]T_N(s,k-1)
        +2 (s+1)^2-d^2/N^2 T_N(s+1,k-1) ,
        
        which follows from master and tnsk-1.
        
        
        
         §.§ Generating function
        
        
        In order to calculate the mean value and the variance of the reaction time t=k/N 
        we introduce the generating function
        _N(s,z)=∑_k=1^∞ z^kT_N(s,k) .
        
        
        The evolution from s'=n to s'=s in figure <ref>(a) proceeds through a succession of steps
        where the system remains for some time in a state 
        with s'_A+s'_B=2s' until two particles annihilate and s'→ s'-1.
        We associate with such a step the generating function for its lifetime corresponding 
        to the diagrams of figure <ref>(b):
        _N(s',z)   ={1+z[1-2s'^2-d^2/N^2]+⋯
        +z^l[1-2s'^2-d^2/N^2]^l_l updates without annihilation
        +⋯}2zs'^2-d^2/N^2_annihilation
           =2z(s'^2-d^2)/N^2-z[N^2-2(s'^2-d^2)] .
        
        The generating function for the reaction time, measured in the number of updates, is obtained 
        as the product:
        _N(s,z)=∏_s'=s+1^n_N(s',z)=(n-d)!(n+d)!/(s-d)!(s+d)! (2z)^n-s/∏_s'=s+1^n{N^2-z[N^2-2(s'^2-d^2)]} .
        
        It is easy to verify on this expression that _N(s,1)=∑_k=1^∞ T_N(s,k)=1
        so that T_N(s,k) is properly normalized.
        
        
        
         §.§ Mean value and variance 
        
        The mean value of the reaction time t=k/N is given by:
        t_N(s)=1/N∑_k=1^∞ k T_N(s,k)
        =1/N.∂_N/∂ z|_z=1
        =n∑_s'=s+1^n1/s'^2-d^2 .
        
        When d=0 one obtains
        t_N(s)=n∑_s'=s+1^n1/s'^2=n(H_n^(2)-H_s^(2)) ,
        
        where H_l^(m)=∑_j=1^l1/j^m is a generalized harmonic number.
        When d>0 one may write
        t_N(s)=n/2d∑_s'=s+1^n(1/s'-d-1/s'+d)
        =n/2d(H_n-d-H_s-d-H_n+d+H_s+d) ,
        
        where H_l=H_l^(1) is a harmonic number.
        
        A second derivative gives the mean square value of the reaction time:
        t_N^2(s)   =1/N^2∑_k=1^∞ k^2 T_N(s,k)
        =1/N^2.∂/∂ z(z∂_N/∂ z)|_z=1
           =n^2(∑_s'=s+1^n1/s'^2-d^2)^2
        +n^2∑_s'=s+1^n1/(s'^2-d^2)^2-
        1/2∑_s'=s+1^n1/s'^2-d^2 .
        
        Since the first term in the last expression is t_N(s)^2 the variance is given by:
        Δ t_N^2(s)=
        n^2∑_s'=v+1^n1/(s'^2-d^2)^2-1/2∑_s'=s+1^n1/s'^2-d^2 .
        
        When d=0 one obtains:
        Δ t_N^2(s)=
        n^2∑_s'=s+1^n1/s'^4-1/2∑_s'=s+1^n1/s'^2=n^2(H_n^(4)-H_s^(4))-1/2(H_n^(2)-H_s^(2)) .
        
        When d>0dtn2s-1 can be rewritten as:
        Δ t_N^2(s)   =n^2/4d^2∑_s'=s+1^n
        [1/(s'-d)^2+1/(s'+d)^2]-1/4d(1+n^2/d^2)
        ∑_s'=s+1^n(1/s'-d-1/s'+d)
           =n^2/4d^2(H_n-d^(2)-H_s-d^(2)
        +H_n+d^(2)-H_s+d^(2))-1/4d(1+n^2/d^2)
        (H_n-d-H_s-d-H_n+d+H_s+d) .
        
        
         §.§ Scaling limit
        
        In the scaling limit the probability distribution T_N(s,k) leads to three different 
        probability densities,
        depending on the values of x=s/n and y=d/n. We now study these different cases.
        
        
        
          §.§.§ x>y≥ 0.
        
        
        
        
        We first consider the case where both n,s,d→∞ for fixed values 
        of the ratios x>y≥0. 
        Using in tns-3 the asymptotic expansion for harmonic numbers
        
        H_nα=ln n+lnα+γ+(n^-1) ,
        
        where γ is Euler's constant, gives
        t_N h_y(x)=1/2yln[(1-y)(x+y)/(1+y)(x-y)] ,
        
        which is the mean-field expression tmf-1. It reduces to
        t_N1/x-1 ,
        
        when y→0.
        
        Using the following expansion for generalized harmonic numbers
        
        H_nβ^(2)-H_nα^(2)=∑_j=nα+1^nβ1/j^2
        =∫_nα^nβdj/j^2+(n^-2)
        =1/n(1/α-1/β)+(n^-2)
        
        as well as hn, the scaling limit of the variance follows from dtn2s-3 and reads:
        Δ t_N^2   χ_y(x)/n
        χ_y(x)   =1/2y^2(x/x^2-y^2-1/1-y^2)
        -1/4y(1+1/y^2)
        ln[(1-y)(x+y)/(1+y)(x-y)] .
        
        When y→0 one obtains:
        χ_0(x)=1/3x^3-1/2x+1/6 .
        
        Here too the fluctuations are small and t is a self-averaging variable when x>y.
        
        The dependence on x=ρ_A+ρ_B of the mean value and the variance of the 
        reaction time t is shown in figure <ref> for y=ρ_A-ρ_B=0 and 
        figure <ref> for y=1/8. A good collapse of the finite-size data is obtained. 
        
        
        
        
        
        The scaling of the variance with n suggests the definition of the 
        following scale-invariant and centered time variable
        θ(s,k)=n^1/2(t-t_N) ,     t=k/2n
        where t_N given by tns-4 depends on s through x=s/n.
        Accordingly we define the probability density as:
        (x,θ)2n^1/2T_N(s,k) .
        
        We proceed as in section <ref> [Except that here the expression
        of t_N is already known.] and solve the master 
        equation mastert in the scaling limit. 
        
        
        We replace T_N by ,  
        make a Taylor expansion of  to second order in s and k and rewrite 
        the coefficients and the derivatives in terms of the new scaled variables.
        In this way an expansion in powers of n^-1/2 is obtained. The coefficients of n^0  
        and n^-1/2 vanish identically. The first non-vanishing contribution is coming at order n^-1 
        and leads to the partial differential equation:
        ∂/∂ x=1/2(x^2-y^2)(1/2-1/x^2-y^2)∂^2/∂θ^2 .
        
        Assuming that  depends on x only through the reduced variance χ_y(x)
        given by dtn2-1 the partial differential equation transforms into the 
        diffusion equation
        ∂/∂χ_y=1/2∂^2/∂θ^2 .
        
        Thus the fluctuations are Gaussian:
        (x,θ)=^-θ^2/[2χ_y(x)]/√(2πχ_y(x)) ,     x>y≥0 . 
        
        
        The Gaussian  behaviour of the reaction time is shown for different values of 
        x=ρ_A+ρ_B in figure <ref> for ρ_A=ρ_B and figure <ref>
        for ρ_A-ρ_B=1/8. The finite-size results were obtained by iterating the master 
        equation master, storing the data at each update for a given value of s=nx+1 
        and using tnsk-1.
        
        
        
          §.§.§ x=y>0, s≥ d.
        
        We study now the case where s=d+u with d=ny, s_A=(n) and s_B=u=(1), 
        i.e. when the reaction is close to completion or complete. 
        The mean reaction time which follows from tns-3t_N=1/2y[ln(2ny1-y/1+y)
        +γ-H_u]+(n^-1) ,
        
        has a slow logarithmic growth with n which in the scaling limit yields:
        t_N/ln n=1/2y+[(ln n)]^-1) .
        
        The variance in dtn2s-3 gives
        Δ t_N^2=1/4y^2[ζ(2)-H_u^(2)]
        +(ln n/n) .
        
        The fluctuations of the reaction time t are stronger when the reaction is completed 
        or close to completion. The finite-size data collapse for the mean value 
        and the variance as a function of y is shown in figure <ref>.
        
        In this regime the ratio √(Δ t_N^2)/t_N decreases slowly 
        as (ln n)^-1
        which suggests a new type of statistics. It will be obtained by taking the scaling limit
        directly on the probability distribution T_N(s,k) in tnsk-2.
        Since the variance does not depend on n we define a centered time variable as:
        θ'=2y(t-t_N) ,     t=k/2n .
        
        The factor 2y is suggested by the form of t_N given by tn-2.
        With this definition one obtains:
        
        k=n/y[θ'+ln(2ny1-y/1+y)+γ-H_u] .
        
        Thus one defines the probability density as:
        '(u,θ')n/yT_N(d+u,k) .
        
        With the change of summation variable r=j+d+u+1 tnsk-2 leads to
        n/yT_N(d+u,k)=1/u!∑_j=0^n-d-u-1(-1)^j/j!
        A_jB_jC_jD_j(1-E_j)^k-1 ,
        
        where, in the scaling limit with d=ny:
        
        A_j=d+j+u+1/d=1+j+u+1/ny≃1 ,
        
        B_j=(n-d)!/(n-d-j-u-1)!
        =∏_i=0^j+u(n-d-i)≃(n-d)^j+u+1 ,
        
        C_j=(n+d)!/(n+d+j+u+1)!
        =[∏_i=1^j+u+1(n+d+i)]^-1≃(n+d)^-(j+u+1) ,
        
        D_j=(2d+2u+j+1)!/(2d+u)!
        =∏_i=1^j+u+1(2d+u+i)
        ≃(2d)^j+u+1 ,
        
        E_j=(d+j+u+1)^2-d^2/2n^2
        =2d+j+u+1/2n^2 (j+u+1)
        ≃y/n (j+u+1) .
        
        
        Taking into account the expression of k in k one obtains:
        (1-E_j)^k-1≃(2ny1-y/1+y)^-(j+u+1)^-(j+u+1)(θ'+γ-H_u)
        Thus in the scaling limit nytn gives
        '(u,θ')=1/u!^-(u+1)(θ'+γ-H_u∑_j=0^∞[-^-(θ'+γ-H_u]^j/j! ,
        
        which is the generalized Gumbel distribution <cit.>:
        '(u,θ')=1/u!exp[-(u+1)(θ'+γ-H_u)
        -^-(θ'+γ-H_u)] .
        
        It crosses over to the Gaussian in txtheta-2 when u≫1. 
        The matching of the two probability densities is studied in appendix E.
        The collapse of the finite-size data on the generalized Gumbel distribution is shown 
        in figure <ref> for different values of u and y=ρ_A-ρ_B=1/2.
        The finite-size data were obtained by iterating master and using tnsk-1.
        
        
        
        
          §.§.§ x=y=0, s≥ d.
        
        
        
        
        Finally we consider the case when s=(1)≥ d. According to tns-2 
        and dtn2s-2 in the scaling limit the mean value and the variance of the 
        reaction time behave as 
        t_N/nζ(2)-H_s^(2) ,    Δ t_N^2/n^2ζ(4)-H_s^(4) ,
        
        when d=0 whereas tns-3 and dtn2s-3 lead to
        t_N/n1/2d(H_s+d-H_s-d) ,  Δ t_N^2/n^21/4d^2[2ζ(2)-H_s-d^(2)-H_s+d^(2)+1/d(H_s-d-H_s+d)],
        
        when d>0. The finite-size data collapse is shown in figure <ref> for s_A=s_B 
        and figure <ref> for d=(s_A-s_B)/2=4.
        
        The mean value and the standard deviation are both growing as n, 
        thus the reaction time is a strongly fluctuating random variable 
        when the reaction is almost complete. In the following we use the 
        scale-invariant time variable
        θ”=t/n=k/2n^2 ,
        
        and define the associated probability density as:
        ”(s,θ”)2n^2T_N(s,k) .
        
        
        
        
        In the expression tnsk-2 of T_N(s,k) one may write:
        r(r+s)!/(r-s-1)!=(-1)^s+1∏_m=0^s(m^2-r^2) .
        
        Furthermore, in the scaling limit, one obtains:
        (1-2r^2-d^2/N^2)^k-1≃^-(r^2-d^2)θ” ,    (n-d)!(n+d)!/(n-r)!(n+r)!=∏_j=d^r-1n-j/n+j+1≃1 .
        
        Thus the probability density is given by: 
        ”(s,θ”)=2 ^d^2θ”/(s-d)!(s+d)!∑_r=s+1^∞(-1)^r∏_m=0^s(m^2-r^2) ^-r^2θ” .
        
        
        The collapse of the finite-size data on ”(s,θ”) is shown at different 
        values of s for s_A=s_B in figure <ref> and for d=(s_A-s_B)/2=4 in figure <ref>.
        The finite-size data were obtained by iterating master and using tnsk-1.
        
        The asymptotic behaviour for θ”≫1 (see figure <ref>) is governed 
        by the first term in the sum and reads:
        ”(s,θ”)≃2(s+1)(2s+1)2s s+d^-[(s+1)^2-d^2]θ” ,
            θ”≫1 .
        
        
        In order to study the asymptotic behaviour when θ”≪1 it will be convenient to 
        re-express ”(s,θ”) in terms of Jacobi theta functions. 
        First let us notice that the product in tstheta”-2 vanishes for r≤ s 
        so that the sum can start at r=1 instead of s+1. The product can then be replaced 
        by the differential operator ∏_m=0^s(m^2+d/dθ”) so that: 
        ∑_r=s+1^∞(-1)^r∏_m=0^s(m^2-r^2) ^-r^2θ”
        =∏_m=0^s(m^2+d/dθ”)∑_r=1^∞(-1)^r^-r^2θ” .
        
        Making use of the identity (<cit.> p 463)
        ∑_r=1^∞(-1)^rq^r^2cos(2rz)=ϑ_4(z,q)-1/2 ,
        
        one obtains
        ”(s,θ”)=^d^2θ”/(s-d)!(s+d)!∏_m=0^s(m^2+d/dθ”) ϑ_4(0,^-θ”) ,
        
        where the Jacobi theta function can be re-written as (<cit.> p 475):
        ϑ_4(0,^-θ”)
        =√(π/θ”)ϑ_2(0,^-π^2/θ”)
        =2√(π/θ”)∑_r=0^∞^-(r+1/2)^2π^2/θ” .
        
        Finally the probability density takes the following form:
        ”(s,θ”)=2^d^2θ”/(s-d)!(s+d)!∏_m=0^s(m^2+d/dθ”)√(π/θ”)∑_r=0^∞^-(r+1/2)^2π^2/θ”.
        
        
        The leading contribution when θ”≪1 (see figure <ref>) 
        comes from the derivatives of the exponential in the first term of the sum so that:
        ”(s,θ”)≃2√(2)/(s-d)!(s+d)!(π/2θ”)^2s+5/2^-π^2/(4θ”) ,    θ”≪1 .
        
        
        § CONCLUSION
        
        
        The two-species diffusion-annihilation process A+B→ Ø  has been 
        studied on the fully-connected lattice with size N for either equal or 
        different numbers of particles, s_A and s_B, in the initial state where s_A+s_B=N=2n.
        Exact probability distributions, S_N for s=(s_A+s_B)/2 at a given time and 
        T_N for the reaction time t needed to reach a given number of surviving particles, have been 
        obtained by solving the master equation. The finite-size scaling behaviour of the 
        mean value and the variance of s and t has been determined using a generating 
        function approach. 
        
        In the scaling limit, the fluctuations of the number of particles around their mean values, 
        given exactly by mean-field theory, are weak and Gaussian.
        The statistical properties of the reaction time display 
        three different regimes. When both s_A and s_B are (n) the fluctuations 
        of the reaction time, as measured by the ratio R=√(Δ t_N^2)/t_N,
        decay as n^-1/2. The fluctuations are Gaussian and mean-field theory is exact. 
        For unequal initial densities, in the vicinity of the absorbing state when
        s_A=(n) and s_B=(1), one obtains R=[(ln n)^-1] thus the fluctuations 
        are “marginally weak”. They are governed by a generalized Gumbel distribution, 
        indexed by s_B, which crosses over to the Gaussian density with increasing values of s_B.
        For equal or almost equal initial numbers of particles and when the system is
        close to the absorbing state, i.e. when s_A=(1)≥ s_B≥0, one obtains R=(1).
        The reaction time is then strongly fluctuating. Its probability density involves 
        an alternating infinite series which can be considered as resulting from the applications 
        of a product of s+1 first-order differential operators to a Jacobi theta function.
        
        A generalized Gumbel distribution has been recently shown to govern the fluctuations of the 
        covering time of the fully-connected lattice, i.e. the time needed for a random 
        walker to visit almost each site at least once <cit.> (see also <cit.>). 
        Although the two problems present some analogies (they share the same value of R 
        and the same probability density), the mean values and the variances scale differently 
        with N, the number of random walkers is constant for the covering time and 
        decreasing for the reaction time.
        
        The extreme value statistics obtained for the reaction time near the absorbing state, 
        for almost equal initial numbers of A and B particles, is quite similar 
        to what was obtained in <cit.> for the coagulation process A+A→ A. 
        Probability densities with nearly the 
        same form of asymptotics govern the behaviour of the squared width of an interface generated 
        by a periodic Brownian motion <cit.> as well as the maximum height of the 1D Edwards-Wilkinson 
        model <cit.>.
        
        It seems reasonable to conjecture that similar extreme value statistics should govern 
        the reaction time when the system is close to its absorbing state for D>D_c=2
        when s_A-s_B=(n) in the initial state and for D>D_seg=4 when s_A-s_B=(1)
        in the initial state.
        
        
        
        
        § EVALUATION OF V_R^(R)
        
        According to sns0, in the initial state:
        
        S_N(n,0)=1=v_n^(n) .
        
        The recursion relation rec gives
        
        S_N(n-1,0)=0=v_n-1^(n-1)+v_n-1^(n)=v_n-1^(n-1)-n^2-d^2/n^2-(n-1)^2 ,
        
        so that:
        
        v_n-1^(n-1)=n^2-d^2/n^2-(n-1)^2 .
        
        In the same way
        
        S_N(n-2,0)=0=v_n-2^(n-2)+v_n-2^(n-1)+v_n-2^(n)
        where
         v_n-2^(n-1)   =-(n-1)^2-d^2/(n-1)^2-(n-2)^2v_n-1^(n-1)
        =-(n^2-d^2)[(n-1)^2-d^2]/[n^2-(n-1)^2][(n-1)^2-(n-2)^2] ,
         v_n-2^(n)   =-(n-1)^2-d^2/n^2-(n-2)^2v_n-1^(n)
        =(n^2-d^2)[(n-1)^2-d^2]/[n^2-(n-1)^2][n^2-(n-2)^2] ,
        
        so that:
        
        v_n-2^(n-2)=(n^2-d^2)[(n-1)^2-d^2]/[n^2-(n-2)^2][(n-1)^2-(n-2)^2] .
        
        After lengthy but straitforward calculations the same procedure leads to:
        
        v_n-3^(n-3)=(n^2-d^2)[(n-1)^2-d^2][(n-2)^2-d^2]/[n^2-(n-3)^2]
        [(n-1)^2-(n-3)^2][(n-2)^2-(n-3)^2] .
        
        These results suggest the following conjecture
        
        v_r^(r)=∏_j=0^n-r-1(n-j)^2-d^2/(n-j)^2-r^2 ,
             r=d,…,n-1 ,     v_n^(n)=1 , 
        
        from which vrr can be deduced.
        
        
        
        
        
        
        § CALCULATION OF Ω_R,0(W) AND ITS FIRST DERIVATIVE AT W=1
        
        According to ordw when d=0Ω_r,0(w)=∑_s=0^r(-1)^r-sr+s sr s2sw^s/r+s .
        
        Using the following combinatorial identity for 
        Legendre polynomials <cit.>
        P_r(2w-1)=∑_s=0^r(-1)^r-sr+s sr sw^s
        or0w-1 can be rewritten as:
        Ω_r,0(w)=2w^-r∫_0^wdP_r(2u-1)/duu^rdu .
        
        An integration by parts gives:
        Ω_r,0(w)=2P_r(2w-1)-2rw^-r∫_0^wP_r(2u-1)u^r-1du .
        
        
        When w=1 one has:
        Ω_r,0(1)=2-2r∫_0^1P_r(2u-1) u^r-1du
        =2-r∫_-1^1(x+1/2)^r-1P_r(x)dx .
        
        According to the identity <cit.>∫_-1^1(1+x)^σ P_ν(x)dx=
        2^σ+1[Γ(σ+1)]^2/Γ(σ+ν+2)Γ(σ-ν+1) ,
        
        the integral in or01-2 vanishes and or01-1 is obtained.
        
        The first derivative of or0w-3 at w=1 gives:
        .dΩ_r,0/dw|_w=1
        =2.dP_r(2w-1)/dw|_w=1-2r
        
        The generating function for Legendre polynomials
        ∑_l=0^∞ P_l(x)u^l=1/√(1-2ux+u^2)
        leads to:
        ∑_l=0^∞.dP_l(x)/dx|_x=1u^l=u/(1-u)^3
        =∑_j=0^∞j+2 2u^j+1 .
        
        Identifying the coefficients of u^r, one obtains
        .dP_r(2w-1)/dw|_w=1=2 r+1 2=r(r+1) ,
        
        and dor01-2 leads to dor01-1.
        
        
        
        
        
        
        § SCALING LIMIT OF S_N/N AND 
        S_N^2/N^2 WHEN D=0
        
        
        According to snk one may write: 
        s_N(k)/n=∑_r=1^n f(r) ,     f(r)
        =2r/n-r∏_j=1^rn-j/n+j_p(r)  (1-2 r^2/N^2)^k_e_k(r) .
        
        Making use of the following expansions in powers of n^-12r/n-r≃2r/n(1+r/n+r^2/n^2) ,
        ln p(r)   =∑_j=1^r[ln(1-j/n)-ln(1+j/n)]
        ≃-2∑_j=1^r(j/n+j^3/3n^3)
           ≃-r(r+1)/n-r^2(r+1)^2/6n^3 ,
        ln e_k(r)=kln(1-2r^2/N^2)
        ≃-2k(r^2/N^2+r^4/N^4)
        ≃-t(r^2/n+r^4/4n^3),
        
        one obtains:
        
        f(r)=2r/n(1+r/n+r^2/n^2)
        exp[-(t+1)r^2/n-r/n-r^2(r+1)^2/6n^3
        -tr^4/4n^3] .
        
        In the scaling limit the Euler-Maclaurin summation formula gives
        ∑_r=1^n f(r)≃∫_0^∞ f(r) dr+1/2 [f(∞)-f(0)]
        +1/12 [f'(∞)-f'(0)]+… ,
        
        with
        
        f(0)=f(∞)=f'(∞)=0 ,     f'(0)=2/n .
        
        With the change of variable u=r^2/nfr and eml lead to:
        ∑_r=1^n f(r)   ≃∫_0^∞du(1+√(u/n)+u/n)
        (1-√(u/n)-u^2/6n
        -tu^2/4n+u/2n)^-(t+1)u-1/6n
           ≃∫_0^∞du[1+u/2n-
        (3t+2)u^2/12n]^-(t+1)u-1/6n .
        
        Using ∫_0^∞du u^a^-bu=a!/b^a+1 one obtains 
        s_N(k)/n as given in snt.
        
        The expression of s_N^2(k)/n^2 differs from snkn only 
        through a factor u=r^2/n in the sum over r. Furthermore the correction of order n^-1
        to the integral in eml now vanishes. Thus we can directly modify the integral in sfr
        to write:
        s_N^2(k)/n^2≃∫_0^∞du[u+u^2/2n
        -(3t+2)u^3/12n]^-(t+1)u .
        
        Finally, the integration leads to the expression given in snt.
        
        
        
        
        
        
        § SOLUTION OF EQUATION DKAPDT
        
        
        With κ_y(t)=α_y(t)β_y(t)dkapdt leads to 
        β_y(α_y'+4g_yα_y)+β_y'α_y+f_y=0 ,     
        f_y(t)=dg_y/dt+1/2(dg_y/dt)^2 ,
        
        which transforms into a system of two first-order differential equations:
        {[ α_y'+4g_yα_y         =0 ,;  β_y'α_y+f_y         =0 . ].
        
        Since g_y(t) in xmean can be written as
        
        g_y(t)=y+G'/G ,     G=1+y-(1-y) ^-2yt ,
        
        the first equation in ab gives:
        α_y(t)=C_1 ^-4∫^tg_y dt'=C_1 ^-4yt/G^4 .
        
        For the second equation one has
        β_y(t)=-∫^tf_y/α_y dt'+C_2
        
        Where according to abf and gt
        f_y(t)=-4y^2(1-y^2) ^-2yt/G^2+8y^4(1-y^2)^2 ^-4yt/G^4 .
        
        The integration in b-1 is straightforward and gives:
        β_y(t)=2y(1-y^2)/C_1[(1+y)^2 ^2yt-(1-y)^2 ^-2yt-4y(1-y^4)t]+C_2 .
        
        The product of a and b-2 leads to:
        κ_y(t)=2y(1-y^2) ^-4yt/G^4[(1+y)^2 ^2yt-(1-y)^2 ^-2yt-4y(1-y^4)t]+C ^-4yt/G^4 .
        
        The initial condition κ_y(0)=0 is satisfied when C=C_1C_2=-8y^2(1-y^2) 
        and kappat-1 is finally obtained.
        
        
        
        § MATCHING OF '(U,Θ') WHEN U≫1 
        WITH (X,Θ) WHEN X≪1.
        
        
        Let us rewrite '(u,θ') in tutheta'-3 as:
        '(u,θ')=^F_a(b) ,   F_a(b)=-ab-^-b-lnΓ(a) ,   a=u+1 ,
           b=θ'+γ-H_u .
        F_a(b) has a maximum at b_0=-ln a where 
         
        F_a(b_0)=aln a-a-lnΓ(a) ,    
        F_a”(b_0)=-a . 
        
        Since a=u+1≫1 one has
        lnΓ(a)   ≃ aln a-a+1/2ln(2π/a) ,    
        H_u≃ln u+γ ,
        
        b-b_0   =b+ln a=θ'+γ-H_u+ln(u+1)≃θ' ,
        
        and the expansion of F_a(b) around b_0 gives:
        
        F_a(b)≃-1/2ln(2π/u)-u/2θ'^2 .
        
        Thus when u≫1 the maximum is amplified, which justifies the approximation, and one obtains:
        '(u,θ')≃^-uθ'^2/2/√(2π/u) ,     u≫1 .
        
        According to theta and theta' one has
        θ'=2y/n^1/2θ ,
        
        and the change of variables yields
        (x,θ)=2y/n^1/2'(u,θ')=
        ^-(4uy^2/n)θ^2/2/√(2π/(4uy^2/n)) ,     x=y+u/n ,
        
        i.e. a Gaussian with variance Δθ^2=n/(4uy^2).
        This expression has to be compared to txtheta-2 when x-y=u/n≪1. It is easy to verify
        that, in this limit, the variance Δθ^2=χ_y(x) in dtn2-1 
        is governed by the first term on the right and reads
        χ_y(x)≃n/4uy^2 ,
        
        as expected.
        
        
        
        § MATCHING BETWEEN JACOBIAN AND GAUSSIAN REGIMES
        
        We were not able to put in evidence this matching at the level of the probability densities. 
        Thus we shall compare the mean values and the variances in the appropriate limits: 
        x≪1 with x>y≥0 for Gauss and s≫1 with s≫ d≥0 for Jacobi.
        
        For the Gaussian density x>y-1 and x>y-2 yields
        t_N≃1/x ,   y=0 ;    t_N≃1/2yln(x+y/x-y) ,   y>0 .
        
        for the mean value and 
        Δ t_N^2≃1/3nx^3 ,   y=0 ;  Δ t_N^2≃1/4ny^2[1/x-y+1/x+y+1/yln(x-y/x+y)] ,   y>0 ,
        
        for the variance.
        
        In the Jacobian regime using the expansions 
        H_s^(l)≃ζ(l)-s^1-l/(l-1) for l>1 and 
        H_s± d≃ln(s± d) in x=y=0-1 gives
        t_N≃n/s ,   d=0 ;    t_N≃n/2dln(s+d/s-d) ,   d>0 ,
        
        for the mean value, in agreement with x>y-4. 
        Finally, given the above expansions, x=y=0-2 leads to
        Δ t_N^2≃n^2/3s^3 ,   d=0 ;    Δ t_N^2≃n^2/4d^2[1/s-d+1/s+d+1/dln(s-d/s+d)]
         ,   d>0 ,
        
        for the variance, in agreement with x>y-5.
        
        
        
        
        § REFERENCES
        
        
        99alcaraz94 Alcaraz F, Droz M, Henkel M and Rittenberg V 1994  Ann. Phys. 230
          250
        
        hinrichsen00 Hinrichsen H 2000  Adv. Phys. 49 815
        
        benavraham00 ben-Avraham D and Havlin S 2000  Diffusion and Reactions 
        in Fractals and Disordered Systems (Cambridge: Cambridge University Press)
        schutz01 Schütz G 2001 Exactly solvable models for many-body systems far from
          equilibrium  Phase Transitions and Critical Phenomena vol 19 ed Domb C
          and Lebowitz J (London: Academic Press) p 1
        
        odor04Ódor G 2004  Rev. Mod. Phys. 76 663
        
        henkel08 Henkel M, Hinrichsen H and Lübeck S 2008  Non-equilibrium phase
          transitions: absorbing phase transitions vol 1 (Heidelberg: Springer)
        
        odor08Ódor G 2008  Universality in non-equilibrium lattice systems
          (Singapour: World Scientific)
        
        krapivsky10a Krapivsky P L, Redner S and Ben-Naim E 2010  A Kinetic View of Statistical
          Physics (New York: Cambridge University Press) p 414
        
        tauber17 Täuber U C 2017  Ann. Rev. Cond. Matter Phys. 8 1
          
        smoluchowski16 Smoluchowski M 1916  Physik. Z. 17 557
        
        burlatsky87 Burlatsky S F and Ovchinnikov A A 1987  Sov. Phys. JETP 65 908
        
        ovchinnikov78 Ovchinnikov A A and Zeldovich Ya B 1978  Chem. Phys. 28 215
        
        toussaint83 Toussaint D and Wilczek F 1983  J. Chem. Phys. 78 2642
        
        kang84 Kang K and Redner S 1984  Phys. Rev. Lett. 52 955
        
        galfi88 Gálfi L and Rácz Z 1988  Phys. Rev. A 38 3151
        
        krapivsky10b Krapivsky P L, Redner S and Ben-Naim E 2010  A Kinetic View 
        of Statistical Physics (New York: Cambridge University Press) p 435
        
        bramson88 Bramson M and Lebowitz J L 1988  Phys. Rev. Lett. 61 2397
        Erratum 1989  Phys. Rev. Lett. 62 694
        
        bramson91 Bramson M and Lebowitz J L 1991  J. Stat. Phys. 65 941
        
        lee95 Lee B P and Cardy J 1995  J. Stat. Phys. 80 971
        
        tauber05 Täuber U C, Howard M and Vollmayr-Lee B P 2005 
         J. Phys.A: Math. Gen. 38 R79
        
        leyvraz92 Leyvraz F and Redner S 1992  Phys. Rev. A 46 3132
        
        benavraham86 Ben-Avraham D and Redner S 1986  Phys. Rev. A 34 501
        
        hilhorst04 Hilhorst H J, Deloubrière O, Washenberger M J and Täuber U C 
        2004  J. Phys.A: Math. Gen. 37 7063
        
        lee84 Lee K and Weinberg E J 1984  Nucl. Phys. B 246 354
        
        ovchinnikov89 Ovchinnikov A A, Timashev S F and Belyi A A 1989  Kinetics 
        of Diffusion Controlled Chemical Processes (Hauppauge: Nova Science)
        
        savara10 Savara A and Weitz E 2010  J. Phys. Chem. C 114 20621
        
        vardeny80 Vardeny Z, O'Connor P, Ray S and Tauc J 1980  Phys. Rev. Lett. 44 1267
        
        turban18 Turban L and Fortin J-Y 2018  J. Phys. A: Math. Theor. 51 145001
        
        ojo2001 Ojo M O 2001  Kragujevac J. Math. 23 101
        
        pinheiro2016 Pinheiro E C and Ferrari S L P 2016  J. Stat. Comp. Sim. 86 2241
        
        riordan79 Riordan J 1979  Combinatorial identities (Huntington, New York: 
        Robert E. Krieger Publishing Company) p 66
        
        gradshteyn80 Gradshteyn I S and Ryzhik I M 1980  Tables of Integrals, 
        Series, and Products (New York: Academic Press) (7.127) p 797
        
        whittaker27 Whittaker E T and Watson G N 1927  A Course of Modern 
        Analysis (Cambridge: Cambridge University Press) p 462 
        
        turban15 Turban L 2015  J. Phys. A: Math. Theor. 48 445001
        
        chupeau15 Chupeau M, Bénichou O and Voituriez R 2015  Nat. Phys. 11 844
        
        foltin94 Foltin G, Oerding K, Rácz Z, Workman R and Zia R 1994  Phys.
        Rev. E 50 R639
        
        edwards82 Edwards S F and Wilkinson D R 1982  Proc. Roy. Soc. London Ser. A 381 17
        
        majumdar04
        Majumdar S N and Comtet A 2004  Phys. Rev. Lett. 92 225501
        
        majumdar05a
        Majumdar S N and Comtet A 2005  J. Stat. Phys. 119 777
        
        majumdar05b
        Majumdar S N 2005  Current Science 89 2076
        
        fortin15 Fortin J-Y and Clusel M 2015  J. Phys. A: Math. Theor. 48 183001
