Non-universal non-equilibrium critical dynamics with disorder

We investigate finite size scaling aspects of disorder reaction-diffusion processes in one dimension utilizing both numerical and analytical approaches. The former averages the spectrum gap of the associated evolution operators by doubling their degrees of freedom, while the latter uses various techniques to map the equations of motion to a first passage time process. Both approaches are consistent with nonuniversal dynamic exponents, and with stretched exponential scaling forms for particular disorder realizations.

We investigate finite size scaling aspects of disorder reaction-diffusion processes in one dimension utilizing both numerical and analytical approaches. The former averages the spectrum gap of the associated evolution operators by doubling their degrees of freedom, while the latter uses various techniques to map the equations of motion to a first passage time process. Both approaches are consistent with nonuniversal dynamic exponents, and with stretched exponential scaling forms for particular disorder realizations.

I. INTRODUCTION
As the concept of dynamic ensemble distribution remains elusive, a way forward to the study of nonequilibrium systems relies partly on prototypical stochastic models [1]. The identification of scaling regimes for these latter at both large time and length scales has permitted to characterize a vast amount of nonequilibrium processes in terms of universality classes [1]. However, their asymptotic dynamics can change completely in the presence of quenched disorder or under space-dependent external forces [2]. This is because particle motion may well remain localized around strong spatial heterogeneities which prevail over stochastic fluctuations [2,3].
In this context, as well as for their value in modelling real non-equilibrium behaviour, reaction-diffusion models have been much studied in recent years, with such methods as real space renormalization, yielding very detailed analytical results at large times [4]. In the thermodynamic limit, some universal and non-universal types of behavior have been reported for these systems [5] and related models [6], yet comparatively little is known on their relaxation times in finite disordered samples. In this work we focus on the dynamic exponents associated to these time scales as well as on finite size scaling aspects which ultimately will turn out to be nonuniversal. Although these processes are here schematized as one dimensional annihilating random walks in a Brownian potential [4,5], we expect them to be relevant in the description of such real cases as exciton dynamics on long disordered polymers, given the experimental success of their homogeneous counterparts [7].
Typically, the microscopic dynamical rules of these models involve hard-core particles which hop randomly and annihilate in pairs on adjacent locations. To determine the effect of quenched disorder on dynamic exponents, here we allow for both varying diffusion and reaction rates in a linear chain. So we consider the following one-step processes. A particle at site l (l + 1) , chosen from L locations, hops with rate h l (h ′ l ) to site l + 1 (l ) provided it is vacant. Independently, reacting attempts take place with rate R l whenever two particles occupy sites l and l + 1. As we shall see in Secs. III and IV, the nonuniversal situation alluded to above will be associated to random orientation changes in diffusion biases h − h ′ . This scenario is illustrated in Fig. 1, where an evolution snapshot of these stochastic rules is depicted for the case of instantaneous reactions under a binary distribution of biases.
Despite this simplicity, such modeling entails a variety of interesting aspects which play a crucial role in the subsequent analysis. These are: the nonequilibrium character of the process, formalisable in terms of a non-Hermitian quantum Hamiltonian [8]; criticality -associated with large (domain) scales, and small energies; disorder -the distributions are the source of the remarkable non universality; particle nonconservation -bringing new features to the non-universal critical phenomena through the interplay between particle and hole sectors.
This work is organized as follows. In Sec. II we present a calculation approach that exploits the quantum analogy and focuses on the spectrum gap (real part) of the associated non-Hermitian Hamiltonian. For numerical convenience the difficulties brought in by the particle nonconservation are remedied by duplicating the degrees of freedom of the original process. Sec. III outlines an analytic treatment (cf. [9,10]), developed from equations of motion mappings combined with renormalisation and transfer matrix techniques generalised to the non-conserving non-equilibrium nature of the system. Among the main features and procedures (briefly presented below) are: states appearing at quantised values of a complex phase obtained from phase steps for each bond and sign alternations at limiting phases; identification and aggregation of phase steps, as generalised random walks, in general biased, in the region between the limiting phases; this makes the determination of phases and their distributions a first passage problem [11] in which both the character (power law or stretched exponential) of the non-universal critical dynamics and in particular the critical exponents are determined by the effective diffusion constant and bias, inheriting their parameter-dependences (non-universality). Sec. IV implements numerically the ideas developed in Sec. II to check out the theoretical results obtained in Sec. III. The computation of finite size scaling properties of the gap yield the dynamic exponents z, thus providing the fundamental scaling between length and time. In addition, numerical simulations are carried out to corroborate the nonuniversal findings for z. Sec. V contains a summarizing discussion along with some remarks on extensions of this work.

II. FERMION DOUBLING
The common starting point for both theoretical and numerical approaches (Secs. III and IV respectively), is the representation of the evolution operator of the stochastic process, described above, as a quantum spin-1 2 non Hermitian Hamiltonian [8], which henceforth we recast as a problem of spinless fermions via a Jordan-Wigner transformation (JW) [12]. Furthermore, these particles become actually free so long as the constraint is imposed for all lattice bonds [8,13]. Specifically, it can be readily checked that in terms of fermion operators C l , C † l , the evolution operator under these conditions reads . According to the parity e iπNC ≡ e iπ P j C † j Cj of the subspace considered, here is a boundary term stemming from the JW transformation under periodic boundary conditions (PBC). The bilinear form imposed throughout by free fermion constraints is already very suitable for theoretical analysis (as given in Sec. III), however the pairing terms of H C would pose severe size restrictions in its numerical analysis, demanding the consideration of matrices growing exponentially with the length of the system. In this respect, notice that a standard Bogoliubov transformation [14] is not effective within our disordered non Hermitian context. Rather, here we follow the strategy of Ref. [15] in formally related systems and proceed to double the degrees of freedom of the original problem by introducing a replica H D of the Hamiltonian H C . The idea is to find a simple (i.e. disorder independent) unitary transformation so as to cancel out all non-conserving particle terms of the replicated system H C+D = H C + H D . To this end, we put forward the following fermion operators with a's and b's coefficients such that all pairing terms are ultimately canceled in the f − g representation. In fact, this can be thought of as a local Bogoliubov transformation in the real space of the double system. After some simple but lengthy algebra, we find that these coefficients should satisfy Among the various solutions we choose, which leave us with a Hamiltonian of real parameters given by As before, here B C+D parallels the boundary terms referred to in Eq. (3), which in this language take the form Therefore, the number of f and g particles is conserved as all pairing terms actually cancel out, provided we restrict attention to subspaces having (−1) NC = (−1) ND . In passing, it is instructive to check how the exact solution of the homogeneous situation is recovered by Fourier transforming the f l and g l operators respectively to a set of new wave fermions f q , g q . To account for the parity conserving terms implicated in B C+D , the q-momenta must be those of the set Q + = {±π/L, ±3π/L, · · · , ±(L − 1)π/L} for the even subspace, while for the odd one the q's should belong to Q − = {±2π/L, ±4π/L, · · · , ±(L − 2)π/L, 0, π}, say for L even. After straightforward algebraic steps, it is readily found that the ordered Hamiltonian H o [ i.e. R l = R, h l = h, and h ′ l = h ′ ∀ l in Eq. (7) ] reduces to Using a simple non unitary 2 × 2 similarity transformation to new F ± q , G ± q fermion like operators [16], the double system is finally diagonalized as Thus, in doubling the degrees of freedom the Λ q excitations of the original problem are restored [13], but now they appear forked in two symmetrical branches Λ q , −Λ −q around a Fock vacuum with "energy" R L rather than with vanishing value, as it would correspond to the stationary state of the original stochastic H C . Hence, in general the doubling transformation (4) should not be regarded as a similarity one, although it preserves all the anticommutation relations which is enough to maintain the spacing of eigenlevels, and in particular the spectrum gap. In the homogeneous case, this scales as 1/L 2 , so the usual dynamic exponent z = 2 is restituted.
For the more interesting disordered situation, if we think of {f 1 , g 1 , · · · , f L , g L } as the 2L single particle states of the replicated Hamiltonian H C+D , then the single excitations of Eq. (7) can be obtained by diagonalizing the real non symmetric five-diagonal matrix (with boundary conditions) where R = j R j . The numerical diagonalization of M ± [17] reveals the emergence, as before, of two symmetrical branches {R + Λ}, {R − Λ} of single excitations around a ground state energy R. Denoting by Λ ± 0 the lowest level excitation in either subspace ( Re Λ ± 0 = min {Re Λ ± } ), the spectrum gap g of the original evolution operator (2) with an initially even (odd) number of reacting-diffusing particles, is finally obtained as g = 2 Re Λ + 0 , with a doubly degenerate Λ + 0 (g = Re Λ − 0 , with Λ − 0 being non degenerate). We will recapitulate these ideas in the numerical diagonalizations of Sec. IV so as to account for the theoretical findings of the next section.

III. PHASE STEPS AS RANDOM WALKS
In what follows we reconsider the non-conserving fermion Hamiltonian (2), but without fermion doubling. One of the eigenvalue equations associated to that operator involves only "hole" state amplitudes a l at each site l. These provide an inhomogeneous term in the other equation, for the "particle" states a + l . Using two-component vectors, the eigenvalue equation involves 2 × 2 matrix coefficients which depend on the excitation "energy" Λ, and on bonddependent biased hopping and bond-independent pair annihilation rates h l , h ′ l , R l , constrained as in Eq. (1), to keep the free-fermion character. The product of the associated 4 × 4 transfer matrices, L l=1 T l , provides the spectrum through its trace, or appropriate matrix elements (depending on boundary conditions), effectively corresponding to quantisation of an aggregated phase. One 2 × 2 part of T l is the transfer matrix t l for the self-contained subspace ("hole sector") provided by the states a l .
From the coupled eigenvalue equations one may equivalently proceed via the corresponding 4 × 4 transfer matrix T l , or via renormalisation using a decimation method, or by mapping ratios of successive amplitudes. Here we mainly employ the last approach. The transfer matrix and decimation approaches are very instructive, in clarifying respectively the role of accumulating complex phases (which provide the spectrum) and the renormalization of quenched probability distributions of such aggregated variables, but in what follows it is sufficient and convenient to begin with the transfer matrix and then move over to the mapping approach.
For orientation and to obtain certain characteristic parameters it is useful to separate the rates (and consequently the transfer matrix itself) into reference uniform, l-independent, "effective medium" parts (h, h ′ = R − h, indicated by the absence of a site label), together with the further random part involving a disorder strength δ, and (binary) random variables ζ l at each bond. Then T l becomes T l ≡ T +∆ l . For any finite size L, it is possible to start from T and build in the disorder terms∆ l completely. For that, it is convenient to work in the representation diagonalising T . Two (i = 1, 2 say) of the four labels on its eigenvalues λ i relate to the "hole" sector part t l of the transfer matrix. The "particle" sector (i = 3, 4) couples to the other sector through the inhomogeneous term referred to above. The eigenvalues are conveniently parametrised as λ 1,2 = h ′ /h exp(±iq), where q satisfies Λ = 1 − 2 √ h ′ h cos q, with a similar form (with Λ changed sign) for λ 3 , λ 4 .
The formal expansion in powers of∆ l gives where A n and B n come respectively from the i = 1, 2, and i = 3, 4 terms in the trace, and each involves sums over n disorder variables. For example, where η(mq) is a periodic function involving separations m i of the participating disorder sites and connected with competing aspects of the wavelength and scales set by the local disorder and bias. B n involves a part similar to A n and a second more complex part generated by the inhomogeneous terms.
It is instructive to consider the approximation to the multiple sum in Eq. (13) coming from taking all m i as large. In this case the q-dependence in the functions η(m i q) goes into an overall exponential factor, and the product of these factors comes out of the multiple sum as a complex phase factor exp(±i Lq) also accompanied by further qindependent factors. The summand is then just the product of disorder variables ζ ... and the sum can be shown to be the coefficient of x n in the generating function L l=1 (1 + xζ l ). Taken together with the other factors that leads to L n=0 A n being the exponential of an aggregated complex phase with real and imaginary parts L l=1 φ l and i LRe (q), where For binary random rates the sum of φ l is like a random walk. The case of most interest is when the two possibilities have different signs (orientations) for the bias h l − h ′ l , so steps occur in both directions. The same phase step (though aggregated in various different ways) occurs in the parts of the corresponding approximation for L n=0 B n coming from the inhomogeneous terms; the remaining "homogeneous part" has negligible phase step.
However, in the approximation the location of states then becomes the same as for the pure system -coming just from i LRe (q) since the disorder dependence of the phase is now confined to the real function φ l . So it will be necessary to generalise the result obtained by taking characteristic m's as large, by retaining the registration between local disorder and aggregating phase which is provided by the functions η(m i q). This is most easily done by working with the fully-disordered transfer matrix (or directly from the original equations of motion). That gives, as follows, information on characteristic scales, domain structures, phase accumulations, and the region of validity of the approximation just discussed. For the phases from L n=0 A n (the simplest illustration) only the 2 × 2 subpart t l corresponding to the hole sector of the (fully-disordered) transfer matrix T l is needed.
The crucial aggregating phase variable is ln x l , where x l = (a l−1 /a l ) is the ratio of complex amplitudes a l on successive sites. These x l ratios satisfy From this it follows that x l = g l (h ′ l−1 /h l−1 ), where g l is such that z l ≡ (1 − g l ) −1 satisfies From this (Mobius) map, and its first iteration yielding z l+2 , it can readily be seen that when (needing Λ small) the variable χ l ≡ ln z l "walks" according to χ l+1 ∼ χ l + ln (h l /h ′ l−1 ) while the complex phase step in a l is ln x l = ln [g l (h ′ l−1 /h l−1 )] ∼ ln (h ′ l−1 /h l−1 ) = φ l−1 , as was found in Eq. (14) using the approximation on L n=0 A n . Near z l = z > (or z < ), the mapped z goes to ∞ (or 0, then 1, giving a vanishing amplitude) and changes sign. This procedure thus provides the region (17) of validity of the approximation on L n=0 A n , and generalises it: the behaviour at z = z < , z > is the proper consequence of the functions η(m i q), replacing the sign changes given by exp(±i Lq) in the approximate result.
In the corresponding discussion for the structure and phases in L n=0 B n (coming from "particle" states, with complex amplitudes a + l ), linearity allows a separation into two parts: (i) without, and (ii) with, the mixing to the states in the "hole" subspace just treated. The discussion for the first part ((i)) is similar to that just given, and leads to a walk step ln (a + l−1 /a + l ) ≡ ln c + l with c + l close to 1 if y l ≡ (1 − c + l ) −1 satisfies a condition like (17), with the same limits z < , z > to leading order in Λ. For small Λ, ϕ l ≡ ln y l walks with the same step as χ l = ln z l , and has the same behaviour at the limiting values. For part (ii), the inhomogeneous term coming from the mixing to the states in the other subspace leads to the occurrence again of z l and x l = (a l−1 /a l ), in terms of which, for z < ≪ z l ≪ z > (so Λ is small and x l ∼ (h ′ l−1 /h l−1 )) , the amplitude ratio is a + l /a + 0 = 1 + with γ l = h ′ l /h l , and A 0 = Γ 0 − a 0 /x 1 . The combination of these results from (i),(ii) gives the essence of the complicated result for the approximated L n=0 B n , together with its validity criterion, and again generalises it as necessary to obtain and quantize the imaginary part of the phase, to enumerate the states: the limiting boundary phase conditions provide the sign reversals in amplitude and correspond to (reflecting and absorbing) boundary conditions on the walks described above.
Since the walk steps of χ l and ϕ l contain disorder we need the probability distributions P (χ l ) and P (ϕ l ), which are the same for small Λ. The distributions depend on the boundary conditions, and are most easily treated in a continuum (diffusion) picture. There, l plays the role of time, and the accumulating χ l or ϕ l plays the role of spatial coordinate. From the walk steps described above, within the limits set by Eq. (17), the distributions satisfy a standard biased diffusion equation having diffusion constant D and bias b given by The boundary conditions take the biased problem to first-passage form [11] in which many reflections of the walk occur before absorption is reached, with probability exponentially small in the traversal "length". The boundary conditions set the limits to the excursion (of the walk), making the traversal "length" (long, as required for the continuum picture, in the critical regime of small Λ). Consequently, the first passage "time" scale L ∼ L, is proportional to exp[(|b|/D) ln (Λ 0 /Λ)] , which can be written in the scaling form Λ −1/z with dynamic exponent In the opposite case of large bias, applying for nearly homogeneous systems, the effects of (Re Λ) 1/2 dominate making z close to 2. From the exponential distribution P (χ) ∝ exp[− ( | b |/D) χ] for χ (and similarly ϕ), the energy distribution is appreciable only at χ = ln z < (Λ). That gives the exponential energy distribution It can be seen that the dynamic exponent z is non-universal because the distribution of the random rates allows both signs of the walk steps and is (in general) parameter-dependent. In particular, z can diverge if the bias b becomes zero (for example, at a critical concentration p c in the binary case studied below in Sec. IV). Then a stretched exponential scaling form emerges as follows. For the unbiased case the first passage problem has standard χ 2 ∝ DL relation between the "length" and "time" scales. That makes [ln (Λ 0 /Λ)] 2 ∝ DL, consequently with β = 1/2 and c = c 0 D 1/2 , with c 0 a numerical constant of order one. In this case (where the bias vanishes) the Gaussian distribution for the aggregated phase produces a gap distribution with exponential tails. Consequently, for b → 0 , the ratio of rms gap to mean gap diverges.

IV. NUMERICAL RESULTS
To check out the above expectations, we now investigate numerically the finite size behavior of the gap of the original evolution operator (2) when averaged over independent disorder realizations. To this end, we diagonalize the M ± matrices referred to in Eq. (11) given in Sec. II, and whose dimensions just grow linearly with the system size. For concreteness and comparison with the analytical results of Sec. III, let us consider a binary distribution of bonds with parameters B 1,2 = {R 1,2 , h 1,2 }, and probability P (B) = p δ B,B1 + (1 − p) δ B,B2 . For a given concentration p, we averaged spectra typically over 6000 samples of chains with lengths 50 ≤ L ≤ 1000. For most p the resulting average gaps scale as g ∝ L −z irrespective of parity, and eventually with non-universal dynamic exponents z = z(p), as it was established by Eq. (21) in Sec. III. This situation is exhibited in Fig. 2, where we show a characteristic case of instantaneous reactions under orientation changes in diffusion biases h − h ′ , these latter satisfying the free fermion constraint (1) imposed throughout. On approaching the critical regime defined by the unbiased diffusion equation referred to earlier on in Sec. III, finite size effects become progressively pronounced until at a critical concentration p c , given by the universal stretching factor β = 1/2 involved in (23). Thus, in nearing p c the stochastic dynamics slows down dramatically which, in line with the approach of the preceding section, is reflected in the abrupt increase of dynamic exponents observed in Fig. 3. As we have seen in Eq. (21), for p → p c these latter can be arbitrarily large in the thermodynamic limit, but in practice they are ultimately cutoff by our available system sizes. In this regard, notice that the trend of slopes obtained near the critical region (e.g. lowermost curve of Fig. 2), constitutes a numerical lower bound for these exponents (filled symbols of Fig. 3).  Fig. 2; the main panel shows (part of) the case p = pc (σ/ g ∼ 2), whereas the inset displays a characteristic non-critical situation (p = 0.5, σ/ g ∼ 0.6). Both results have exponential tails (in bold lines). (b) shows a typical well defined bias case ({0.6, 0.2}, p = 0.5), exhibiting instead a narrow distribution closely following a Gaussian (in dashed line, σ/ g ∼ 0.02). As is shown in the inset, here dynamic exponents are diffusive.
Also, in agreement with the theoretical analysis, it is seen that the distribution of spectral gaps becomes broadest at p c (maximum variance σ). This is illustrated in Fig. 4a, after binning gaps in 800 intervals over 60.000 disorder realizations. By contrast, chains with non-fluctuating bias orientations have rather narrow gap distributions, and typically diffusive dynamic exponents (z = 2), regardless of the value of p. This characteristic scenario already anticipated in Sec. III, is exhibited in Fig. 4b.
Simulations.-To allow for an independent examination of nonuniversal aspects, we also carried out numerical simulations in periodic chains of 200, 400 and 800 sites (without duplication). Averages of particle densities were taken over 10 3 histories of 10 3 samples with binary disorder. Starting from homogeneous densities ρ, our results in Fig. 5 clearly indicate an asymptotic decay of the form ρ L (t) ∝ e −t/τL with a relaxation time τ L ∝ L z . As expected, the slope collapse follows closely the non-universal dynamic exponent already identified by direct diagonalization (z ∼ 1.64, in Fig. 2). In turn, the best scaling behavior for well oriented biases is obtained with standard diffusive exponents, which is also in agreement with our previous numerical results (inset of Fig. 4b).
Finally, it is interesting to comment on the regime t ≪ L z . This latter is addressed in Fig. 6 where the density decay below, at, and above p c is exhibited for much larger systems. After averaging 20 samples over 100 histories each, it is observed that in all cases there is an early stage showing the typical t −1/2 decay. However, at p = p c (p > p c ), an incipient crossover to a slower (faster) decay clearly emerges (cf. [5,6]). It would take further evolution decades to  Fig. 4b), displaying instead a typical diffusive exponent. Squares denote data for L = 800, circles 400, and triangles 200. To compare slopes, the data of these latter two sizes were shifted downwards.
estimate the actual form of such regimes, although their nonuniversal characters are already in line with our findings for z(p). The situation for p < p c remains unclear as the existence of a very slow crossover can not be ruled out.

V. CONCLUSIONS
To summarise, after identifying phase steps, and limiting phases and their relationship to generalised biased random walks, in Sec. III we have outlined how the character (power law or stretched exponential) of the non-universal critical dynamics and in particular the critical exponents and their non-universality are determined by the associated diffusion constant and bias (in the case where the fundamental phase step can have either sign).
Among the specific results given above are: in the general biased case the dynamic exponent z is D/|b| and associated distributions of the energy gap are exponential. For vanishing bias z diverges and the usual critical power law behaviour is replaced by a stretched exponential form, possibly signaling the emergence of a glassy dynamic, while the ratio of rms gap to mean gap diverges. All those findings were corroborated numerically using the doubling fermion approach given in Sec. II and Monte Carlo simulations of Sec. IV.
Concerning occurrence of both positive and negative phase steps, that is connected, in a given configuration of randomness, to a succession of domains of alternating bias, whose boundary sites alternate between being traps or repellers for particles. That causes some particles to wander ineffectively for long times before finding a partner for pair annihilation (see Fig. 1), leading to large or divergent dynamic exponents (see Figs. 2 and 3). Among many important things not covered are: the full parameter-dependence of dynamic exponent z from pure value to divergence where the bias vanishes (e.g. in the case of the binary distribution of rates, with concentration parameter p, where z(p) is not monotonic for p in [0 , p c ] or in [p c , 1]); the effect of open boundaries versus PBC (can be very significant in driven non-equilibrium systems) -particularly here regarding the stretching exponent for zero bias -needs further investigation.