Multiscaling in the sequence of areas enclosed by coalescing random walkers

We address the question whether the sequence of areas between coalescing random walkers displays multiscaling and in the process calculate the second moment as well as the two point correlation function exactly. The scaling of higher order correlation functions is estimated numerically, indicating a logarithmic dependence on the system size. Together with the analytical results, this confirms the presence of multiscaling.


Introduction
Gap scaling is found frequently in the context of scale invariance, such as equilibrium statistical mechanics of phase transitions [1], growth phenomena [2], reaction diffusion processes [3] or self organised criticality [4,5]. In the present context, gap scaling [1] means that the nth moment of the observable s, denoted as s n , in leading order scales like t ∆(1+n−τ ) in some large parameter t. The "gap" refers to the constant difference ∆ between two exponents for n and n + 1. In a stochastic process recently introduced as the "totally asymmetric Oslo model" [6,7], the moments of the area under a random walk along an absorbing wall was found to obey gap scaling with τ = 4/3 and ∆ = 3/2.
If the moments scale asymptotically as power laws (possibly with logarithmic corrections) but the exponents do not increase linearly with n, they are said to display multiscaling. The presence of multiscaling is sometimes confused with the absence of scaling altogether and very few examples of simple dynamical processes are known, which can be shown by exact calculation to display multiscaling. It is therefore highly desirable to find a simple stochastic process, such as Brownian motion, for which multiscaling can be derived analytically. Motivated by a recent study of cluster aggregation [8,9], we study the scaling of "sequential moments" or, more accurately, correlation functions of the area size between the trajectories of coalescing random walkers, illustrated in Figure 1. We hope that future research will make the link between our results and the results recently obtained by Munasinghe et al [10] on the scaling of the n-point correlation function of coalescing random walkers.
Recent years have seen a considerable interest in the distribution of the area size under a random walker trajectory. It has been investigated in the context of "directed rice piles" [7,11], as well as in the context of extreme value statistics [12,13]. In both cases, the object of interest is the area between trajectories of coalescing random walkers with initial spacing x 0 , or, equivalently, the area under the trajectory of a random walker, which starts at time t ′ = 0 at distance r(t ′ = 0) = x 0 away from an absorbing wall. In the next section, the model will be introduced in detail.
Addressing problems of extreme value statistics, Majumdar and Comtet [13] constrain the ensemble to trajectories with r(t ′ = 0) = x 0 and, in addition, with r(t ′ = t) = x 0 at a particular "termination time" t ′ = t in the limit x 0 → 0. In this limit, they derive very elegantly the exact distribution function of the area sizes. Due to the lack of additional scales, and, correspondingly, due to the lack of a dimensionless parameter, in the case x 0 → 0 the scaling of all observables is readily derived from dimensional analysis, i.e. the exponents are immediately known and gap-scaling is a necessity for all existing moments. In contrast, we will be concerned with finite x 0 (so that a priori neither exponents nor type of scaling are known) and an ensemble not constrained by the termination time in the form r(t ′ = t) = x 0 . This problem has been addressed earlier [11] using some of the results in [13] for a calculation of the distribution function of the areas in leading order of x 0 .
Unfortunately, only the first moment, which is trivial, is known exactly. Very little is . . starting at t ′ = 0 with initial spacing x 0 and demarcating areas s 1 , s 2 , . . . between them. The dashed lines indicate the "walls" at t ′ = 0 and the terminating time t ′ = t. Whenever random walkers meet, they coalesce and perform the rest of the walk together.
known about higher moments beyond leading order. We employ analytical calculations, based on the formalism introduced in [7] as well as extensive numerics to overcome this problem. In the course, we determine the exact second moment of the distribution, as well as the two point correlation function, which displays anti-correlations, sometimes interpreted as a hallmark of multiscaling [10,14]. To our knowledge, these are the first non-trivial exact results for the area size distribution under a random walker with x 0 > 0.
At least three correlation functions are required to decide about the presence of multiscaling. So far, we have been unable to perform the calculation of the third moment analytically. We have therefore resorted to Monte Carlo techniques, which are, unfortunately, hampered by slow convergence due to rare events. Nevertheless, the numerics strongly suggests multiscaling.
In the next section, we will introduce the model and its observables. The third section contains most of the analytical and numerical results, which are briefly summarised in the last section. Figure 1 shows a lattice realisation (see below) of a sequence of random walkers i = 0, 1, 2, . . . with diffusion constant D 0 and trajectories r i (t ′ ) which start at t ′ = 0 from r i (t ′ = 0) = ix 0 . Whenever two random walkers meet, they coalesce and the resulting, single walker continues its walk. The walkers are indistinguishable and it is therefore irrelevant whether walkers actually run simultaneously or sequentially, in which latter case they are thought to stop as soon as they intersect the trajectory of another walker. This is the picture we will adopt henceforth.

Model
For the sake of definiteness, all trajectories are considered to run from t ′ = 0 to t ′ = t. If two walker coalesce, their trajectories coincide for the remaining time. The area s i between two walkers i − 1 and i with trajectories r i−1 (t ′ ) and r i (t ′ ) respectively, is a random variable, the moments of which will be denoted by s n , where is the expectation value over all configurations.
Given an individual trajectory, it is impossible to tell whether it is produced by merging walkers or not. Producing the trajectories sequentially, the only trajectory relevant to the fate of walker i is the trajectory i − 1. There is no transient (other than producing the very first trajectory), i.e. stationarity is reached immediately. That means, for example, that the expectation value s 1 s 2 is translation invariant, If the random walkers are unbiased, any sequence of areas (s i , s i+1 , . . . , s i+n ) is, by construction, as likely as the sequence (s i+n , s i+n−1 , . . . , s i ). This can be seen by mirroring the entire sequence which, in case of an unbiased walk, leaves the probabilities of individual trajectories unchanged. In the presence of a bias this is no longer case and so the joint probability of a sequence changes under inversion. Thus, in the presence of a bias, one expects, for example, s 2 1 s 2 = s 2 2 s 1 , while the two averages would be equal in the absence of a bias. Of course, the observable itself might be inversion invariant, for example s 1 s 2 = s 2 s 1 regardless of a bias.
As long as the two paths r i and r i−1 have not merged, the variance of the difference r i −r i−1 corresponds to the variance of a random walker with twice the diffusion constant D = 2D 0 , walking along an absorbing wall. However, while the individual trajectories r i might be biased, i.e. in the presence of a drift, the difference does not have such a bias. It is therefore clear that whatever is said about the (single) area size distribution between random walkers with initial spacing x 0 , biased or not, equally applies to the (single) area size distribution under an unbiased random walker with twice the diffusion constant starting at distance r(t = 0) = x 0 from an absorbing wall. However, the notion of a sequence is missing in the latter case.
The main aim of this work is to calculate the correlation functions s 1 s 2 , s 1 s 2 s 3 , and s 1 s 2 s 3 s 4 and determine their asymptotes, that means the exponent γ n and possibly the polynomial P so that lim t→∞ s 1 s 2 . . . s n t γn P(ln t) = c n > 0 (2) suspecting that γ n is not linear in n and/or that P is non-trivial.
As it turns out, s 1 s 2 and indeed s 1 s n+1 can be derived from s 2 , the calculation of which is detailed in the first part of the result section. The other two moments could only be tackled numerically, presented in the second part of the result section.
Most asymptotes will be determined in the limit of large t. It is important to realise the rôle x 0 plays in this context. The problem is fully parametrised by t, D and x 0 . With the area having the same dimension as t √ Dt, in case of vanishing x 0 the scaling of the expectation value of any product of these areas, say s 1 s 2 3 , is fixed to be a corresponding power of t √ Dt, i.e. s 1 s 2 3 ∝ t 9/2 . In fact, there cannot even be any corrections to this behaviour. However, if x 0 = 0, there is a dimensionless quantity Dt/x 2 0 , which can and, indeed, does change the exponents of the asymptotes and gives rise to corrections to scaling [15].

Results
In this section, we will first discuss in detail the analytical results, starting with a brief review of the method which is used to calculate the second moment s 2 in closed form. After deriving its expansion in powers of t, some shortcomings of the perturbative approach of [7] are discussed. Using these results, the correlation function s 1 s n+1 is derived, which includes the original problem s 1 s 2 as a special case. This part finishes with a discussion on the effect of anti-correlations and a brief section on higher moments.
The second part of the result section is concerned with our numerical results. First, the lattice effects are estimated by a comparison of the numerical estimate for s 1 s 2 with the analytical result. After discussing some technical problems, the results for s 1 s 2 s 3 and s 1 s 2 s 3 s 4 are shown, indicating the presence of logarithmic corrections.

Exact calculations
When calculating correlation functions of the form s 1 s 2 , one of the central observations is that they can be obtained by considering "local moments" of the form s 2 in a system with twice the initial spacing, ‡ which will be discussed in detail below. We therefore start the section with the calculation of the second moment and its implications for earlier results, in particular the approximation scheme used in [7]. Knowing s 2 then allows us to calculate s 1 s n+1 .
All analytical calculations are based on the hierarchy of differential equations introduced in [7], with boundary conditions Starting with n = 0, ψ 0 (x, t; x 0 ) is the probability to find at x a fair random walker along an absorbing wall with diffusion constant D = 2D 0 that started at x 0 . For n > 0, ψ n is the local contribution to the n-th moment of the area under the trajectory, see [7]. For n = 0 one finds, consistently, s 0 = 1, for n > 0 one has ‡ We thank Alan Bray for sharing this insight with us.
The formal solution is the hierarchy which is, as it turns out, not easily calculated in closed form. Only the first moment is calculated immediately It has been shown [7,11] that the moments follow gap scaling, with τ = 4/3 and ∆ = 3/2 and C n dimensionless. Using the hierarchy (11) and the definition of the moments (6), the second moment is where T = T (t; x 0 ) = t x 2 0 /D . The integration of (14) would be fairly straight forward, was it not for the factors y and y ′ .
There is no reason to expect convergence problems from the improper integrals, because the Gaussians (and even more so differences of Gaussians) effectively cut them off. We therefore feel confident when changing the integration order, so that the improper integrals are done first and in reverse order. The integration over y can be done immediately and produces only a factor of the form y ′ 4(τ − τ ′ )π. The integral over y ′ is slightly more complicated, as it is of the form dy ′ y ′2 G(y ′ , τ ′ ; 1), and it produces an exponential and an error function, all multiplied by powers of τ ′ . Remarkably, the resulting integrand for the integration over τ ′ and τ is independent of τ , so that the integral over τ ′ can be replaced by the same expression evaluated at τ with a prefactor T − τ . Due to the presence of pre-factors involving various powers of τ , the final integration over τ produces many different terms: To our knowledge this is the first non-trivial moment of the area distribution under a random walker that has been calculated exactly.
It is very instructive to expand the result in powers of T = T (t; with the terms beyond leading order to be considered corrections to scaling [15]. There are two terms without a factor π −1/2 , namely − 1 6 T and − 1 180 , which are present in the same form, i.e. without any further factor, already in the exact expression (15). What makes these terms very different from all others is that they are not an odd-half power of T . In fact, to any order, all terms apart from −(1 + 30T )/180, are odd-half powers of T . In [7] only the odd-half powers had been anticipated, suggesting that no deviation from this form is possible, but that turns out to be incorrect now.
It is instructive to understand why the systematic, perturbative expansion in [7] misses the two terms. It relies on a perturbative expansion of ψ 0 ( √ µy, µτ ), (10), for large µ at fixed, finite y and τ . The scaling in µ is then inherited by ψ 1 ( √ µy, µτ ) through the integral (11) and by s 2 (µt; x 0 ) through (6). This relation gives rise to the scaling of s 2 (t; x 0 ) by parameterising t in convenient units, t = µt 0 . Along the lines of a renormalisation group calculation [16], µ can be considered the "flow parameter" and t 0 = x 2 0 /D as the normalisation point. Using this procedure for all orders in µ produces the prediction of only odd-half powers of µ.
However, when evaluating the perturbative expansion of ψ 0 ( √ µy, µτ ) in powers of µ in the integrals producing ψ 1 and s 2 , one violates the premise that y and τ are finite, as in the integral (11) τ runs from 0 to a finite value and the integral over y is to be evaluated at divergent upper bound. As a result, one finds integrals of the form . .. For i = 0, 1 the integrals are convergent, producing terms of order T 5/2 and T 3/2 respectively, however for i ≥ 2 one (i = 2) or both (i > 2) integrals diverge in the lower limit. Ignoring this divergence, i.e. only keeping the upper limit, produces the correct amplitudes nevertheless. We have tested the integrals for i = 0, 1, . . . , 5 by direct evaluation and comparison to the exact result in the form of the Taylor expansion (17), of course, with the terms of order T 1 and T 0 missing. In general, the integrals for the perturbative expansion give the following odd-half powers: The terms 1/(3/2 − i) and 1/(5/2 − i) are reminiscent of the integration T 0 dτ τ 0 dτ ′ τ ′1/2−i . As long as none of the coefficients or the error term is divergent, the series expansion can be considered a simple Taylor series in µ and therefore in particular the leading order terms are reliable, as in an asymptotic expansion.
Similarly, one can use the perturbative expansion scheme to determine the amplitude of higher orders of s 3 . We can confirm agreement with [11] of the amplitude of the leading order T 4 , expecting the next to leading order T 3 to be correctly predicted in the perturbative approach as well. Following the integration for s 2 , now three integrals over τ might diverge and one might speculate that the first term missed in the perturbative expansion is T 5/2 and the (actually divergent) amplitude for T 2 to be correct, a T 3/2 to be missed, T 1 to be correct, T 1/2 missed again, and terms of order T 0 and lower to be correct again.
We now turn back to the expansion of the exact result, Equation (15). When it comes to the large distance behaviour of the correlation function, the second moment s 2 (t; x 0 ) will be evaluated for small T . Expanding Equation (15) for small T gives Again, the leading, polynomial orders, which together with the pre-factor give just s 2 (t; x 0 ) = x 2 0 t 2 + (2/3)Dt 3 + . . ., are solely due to the two unusual terms with integer powers identified in (15). However, it is consistent with the method in [7] being a large T expansion, that they would have been missed there. All other terms in (19) are exponentially suppressed with a remarkably high leading power in the polynomial. In fact, the first six terms, starting with T 1/2 , all turn out to have vanishing amplitude.
and more generally with S = s 2 + s 3 + . . . + s n−2 for n ≥ 3 (where S = 0 at n = 3) Figure 2. Pruning the network of trajectories, such as the one shown in figure 1, by leaving only every nth trajectory produces an ensemble of walkers with initial spacing nx 0 (trajectories removed are shown dotted). Pruning is identical to merging the enclosed areas, in this example n = 3 so that s ′ 1 = s 1 + s 2 + s 3 , s ′ 2 = s 4 + s 5 + s 6 ,. . . . Such a joining/pruning scheme, mapping an ensemble with initial spacing x 0 onto an ensemble with initial spacing nx 0 is exact, even on the lattice.
for n ≥ 0. Given Equation (15), the second moment u n (t; x 0 ) is easily calculated because u n (t; x 0 ) = u 1 (t; nx 0 ) = s 2 (t; nx 0 ). This is illustrated in Figure 2: The network of trajectories is pruned by joining n consecutive areas s 1 , . . . , s n to one large area. Thereby a sequence of larger areas is constructed, as if produced by random walkers with initial spacing nx 0 . The exact two point correlation function therefore is for n ∈ Z, with the exact expression given in (15). The modulus guarantees the validity even for n ≤ 0, although the above results for s 2 (t; x 0 ) have been derived using x 0 > 0. In the following, the asymptotes of s 1 s n+1 are examined. We first consider the behaviour of s 1 s n+1 (t; x 0 ) in the limit of large t at fixed n. By dimensional analysis, s 2 (t; x 0 ) consists of a pre-factor (x 3 0 /D) 2 and a dimensionless function in T = Dt/x 2 0 . According to the Taylor expansion (17), a term of order T 5/2 gives rise to a contribution in total linear in x 0 , and lower order terms in T correspond to larger powers in x 0 . Under the provision of large T according to Equation (17) the powers are x 0 , x 3 0 , x 4 0 , x 5 0 , x 6 0 , x 7 0 , x 9 0 , . . .. When taking the lattice Laplacian according to (23), Equation (17) is evaluated term by term for x 0 → nx 0 and x 0 → (n ± 1)x 0 , which enters in the pre-factor (x 3 0 /D) 2 as well as in T (t; x 0 ) = tD/x 2 0 . Since (n + 1) − 2n + (n − 1) = 0, the first power of x 0 to contribute therefore is x 3 0 , so that in the limit of fixed n > 0 and large t s 1 s n+1 (t; x 0 ) again in terms of T = T (t; x 0 ) = tD/x 2 0 . § As discussed below, s 1 s n+1 − s 2 < 0 for all T greater than some threshold and n = 0.
Of course, in the limit of large n at fixed t, one expects asymptotic independence, This large n behaviour can be obtained from the small T expansion, Equation (19). The polynomial part T 2 together with the pre-factor (x 3 0 /D) 2 contributes the expected constant (tx 0 ) 2 , while T 3 (x 3 0 /D) 2 is independent of x 0 and its contribution vanishes therefore after operating with the lattice Laplacian. The exponential prefactor exp(−1/(4T (t; nx 0 ))) = exp(−n 2 /(4T (t; x 0 ))) renders the contribution from n−1 exponentially more relevant than from n or n + 1: in the limit of large n at fixed t.

Anti-correlations
The Laplacian structure (22) gives rise to rather peculiar behaviour of the error of the numerical estimate of s . In a naïve implementation one probably takes a series of N samples, s 1 , s 2 , . . . , s N and estimates s by which is indeed an unbiased estimator, s = s [17], independent of N. In the following, we will use an over-bar to indicate numerical estimates, for example s for the estimate of s . Considering its error, for large N one arrives at a picture as cartooned in Figure 3: The walkers start from r 0 (t ′ = 0) = 0, r 1 (0) = x 0 , r 2 (0) = 2x 0 , . . . , r N (0) = Nx 0 , so that for sufficiently large N the very first one, r 0 , is extremely unlikely to interact with the very last one, r N . The total area between r 0 and r N is exactly the sum over all areas, Assuming that the first and the last are in fact independent, suggests that the total area is well modelled by a Langevin type approach r 0,N (t) − r 0,N (t = 0) = t 0 dt ′ η 0,N (t ′ ) § To avoid confusion henceforth the short hand T (without arguments) is used whenever it is to be evaluated at t and x 0 , i.e. where any n dependence has been removed, with two independent noise sources η 0,N (t) with vanishing mean and correlator η 0, . This correlator has variance r 2 0,N (t) = 2Γ 2 t = 2D 0 t, so that Γ 2 = D 0 = D/2. Assuming r 0 and r N are independent, the variance of the numerical estimate is 2Dt 3 /(3N 2 ), which means that the error in the estimate from the correlated sequence s 1 , . . . , s N vanishes faster than in the fully uncorrelated case, where the variance decays like σ 2 (s)/N = 32D 1/2 t 5/2 x 0 /(15N √ π)+. . ., with σ 2 (s) ≡ s 2 − s 2 .
This surprising behaviour is explained by anti-correlations. To see this more clearly, we point out the "sum rule" where the equality of the first and second line is reminiscent of the Laplacian property (22). When considering large N, comparing to (19) confirms, and more importantly  sharply from at n = 0, Equation (17), to at fixed n > 0, see (24). Of course, in the limit of large n asymptotic independence prevails and only exponentially decaying correlations are left 4T see (26). Figure 4 shows s 1 s 1+n together with its numerical estimate.

Other moments
In the pursuit to determine whether the correlation functions s 1 , s 1 s 2 , s 1 s 2 s 3 , . . . display gap scaling or multiscaling in T , any third correlation function, other than s ∝ T , (12), or s 1 s 2 ∝ T 3/2 , (24), is to be determined. Along the same lines as above and using translational invariance, one finds with the leading order of s 3 ∝ T 4 known from [7,11]. In case of an unbiased random walk, the second line simplifies further, for example s 2 1 s 2 = s 1 s 2 2 . Unfortunately, the above expression does not suffice to determine the leading order of s 1 s 2 s 3 : On the right hand side of Equation (35) (s 1 + s 3 ) 3 is unknown, on the right hand side of Equation (36) effectively the same term is unknown, namely s 2 1 s 3 + s 1 s 2 3 . The remaining terms can be determined using the form suggested above, where only the leading fourth order has so far been proven to exist with C 1 = 15/8. In particular, and similarly for the leading order of which is (35) with all unknown terms on the left hand side and all known terms on the right hand side. So, if s 1 s 2 s 3 has leading order higher or equal to (s 1 + s 3 ) 3 , it must be T 4 , given that the leading order always has positive amplitude. If it is of lower order, it might not appear on the right hand side at all, since it might be cancelled by a negative amplitude from (s 1 + s 3 ) 3 . This is potentially a very useful numerical criterion, because if s 1 s 2 s 3 / (s 1 + s 3 ) 3 does not converge to 0, it means that s 1 s 2 s 3 has leading order T 4 ; from numerics, however, it seems that the ratio asymptotically vanishes.

Numerics
All numerical simulations have been done on the lattice, all analytical calculations in the continuum. On the lattice, in every time step the random walker takes one step in the up or down direction with probabilities p and q ≡ 1 − p respectively. The resulting diffusion constant is calculated from the variance after n steps, 4npq = 2D 0 t with t = n∆t and units chosen so that ∆t = 1 and D = 2D 0 = 4pq. In most of our simulations, the walker was unbiased, p = q = 1/2, implying D = 1, and the initial spacing is chosen to be minimal, x 0 = 2. ¶ This choice could, in principle, introduce discretisation effects, which are discussed below.
The random number generator used throughout this study is the "Mersenne Twister" [18]. ¶ Note that x i (t ′ ) is even at even times and odd at odd times.
The central objective of the numerical approach was to determine the large T behaviour of s 1 s 2 s 3 and s 1 s 2 s 3 s 4 . As mentioned earlier, the most naïve implementation estimates s 1 s 2 . . . s n+1 by which, however, introduces correlations. To avoid that, one can produce independent realisations at a significantly higher cost: While N trajectories produce N −n correlated samples, they produce only N/(n + 1) uncorrelated samples, using a new set of n + 1 trajectories for every realisation of s 1 , s 2 , . . . , s n . Obviously, a correlated sample, for instance s 1 s 2 , s 2 s 3 , . . . , s 7 s 8 (N = 9, n = 2), contains multiple uncorrelated sub-samples, for example s 1 s 2 , s 4 s 5 , s 7 s 8 but also s 2 s 3 , s 5 s 6 etc. The correlated sample therefore converges at least as fast as its uncorrelated sub-samples, however, possibly including anti-correlations. It is very important to produce large samples, because of "exponentially rare but important events" [19]. For example, estimating s 1 s 2 s 3 we found an exceptionally large event where s 1 s 2 s 3 was almost 10 12 times bigger than the average s 1 s 2 s 3 up to that point, estimated from about 3.4·10 11 sequential samples, for x 0 = 2 and t = 2 17 . At the end of the simulations, the inclusion of the extreme event still increased the average s 1 s 2 s 3 by a factor of 3.7.
In addition to the statistical error, there are errors from the discretisation and, closely related, the size of x 0 . Ideally, in a numerical simulation of a continuous random walk, x 0 and t were chosen as large as possible, which is, however, computationally very costly. To estimate the influence of lattice effects due to small x 0 , we compared s 1 s 2 (t; x 0 ) to the exact result (15) for different x 0 = 2, 4, 8, 32. The result is shown in Figure 5 in the form s 1 s 2 / s 1 s 2 versus t. As expected, it suggests large x 0 to avoid lattice effects, which clashes with the demand for large T = tD/x 2 0 to determine large time asymptotes in T . Comparing the product of the square of the statistical error and the CPU-time spent on the results for different x 0 indicates that a sequential, i.e. correlated, see (40), simulation of x 0 = 2 is most efficient.
For technical reasons, the error bars shown in the graphs are estimated from the variance of the estimator for small subsamples. For example, the independent samples were created from 1800 individual runs, each consisting of about 7.2 · 10 6 to 61 · 10 6 independent samples totalling to 4.1 · 10 10 samples, produced by fair random walks as described above. However, while s 1 s 2 suggests satisfactory agreement between discrete numerical and continuous analytical results, as mentioned above, rare events made it very hard to estimate s 1 s 2 s 3 or higher correlation functions. The most obvious way to overcome this problem is to use importance sampling [20] which, however, requires a method to produce these rare, important events and a control over their frequency. We decided to introduce a repulsive potential between the walkers which therefore run simultaneously rather than sequentially. Having no analytical reference other than s and s 1 s 2 , it is difficult to assess the quality of this approach. Given the strong impact Sfrag replacements Figure 5. The ratio of the numerical estimate s 1 s 2 and the theoretical value s 1 s 2 . The number of samples varies between 4 · 10 10 and 1.16 · 10 12 . The symbols are shifted relatively to each other to reduce clutter, the error bars span two standard deviations in total. The dashed line marks perfect agreement between numerics and theory. Data marked (i) is from using independent samples, data marked with (s) comes from sequential, correlated runs, Equation (40). The data points marked with arrows have been produced by importance sampling, (r). of rare events on s 1 s 2 s 3 , the quality of the numerical estimate s 1 s 2 compared to s 1 s 2 has only very limited indicative power for the quality of s 1 s 2 s 3 . It is remarkable how large the error of the importance sampling scheme is for s 1 s 2 compared to the naïve methods, see Figure 5, and how small it is for s 1 s 2 s 3 , see Figure 6. Moreover, the numerically estimated error bar of any estimator is expected to suffer from the same problems as the estimator itself, i.e. too small or wrongly biased samples.
The key result for s 1 s 2 s 3 is shown in Figure 6 in the form s 1 s 2 s 3 /[(x 3 0 /D) 3 T 3/2 ], where we divide by (x 3 0 /D) 3 to render the result dimensionless and by T 3/2 to remove the suspected leading order. It remains somewhat inconclusive whether the fit to a logarithm in the intermediate region breaks down for the sequential runs at large T only because of a lack of statistics. The latter is, however, consistent with the observation that underestimation of s 1 s 2 s 3 is more than twice as frequent as overestimation even when averaging over about 4·10 8 samples. The estimates produced from the importance sampling method seem to confirm the presence of the logarithm even for the largest t = 131072.
From the present result, we can quite confidently rule out gap scaling, which would PSfrag replacements Figure 6. The numerical estimate s 1 s 2 s 3 in the form s 1 s 2 s 3 /[(x 3 0 /D) 3 T 3/2 ] in a loglin plot. Symbols corresponding to different initial spacings x 0 are shifted relatively to each other to reduce clutter. Data marked (i) is from independent samples, (s) from sequential samples, (r) from importance sampling with repulsive potential. The straight line is a fit to a logarithm a ln(T /b) to guide the eye. The importance sampling result (see arrows) fully agrees with the logarithmic fit in the intermediate region of the naïve scheme. It remains unclear why for large T the results from independent samples deviate so significantly from the logarithmic behaviour. require s 1 s 2 s 3 /t 2 to converge in the limit of large t to a non-vanishing value. A direct measure for the validity of gap scaling is a moment ratio of the form [6] which converges to a non-vanishing value if the observables obey gap scaling, irrespective of the value of the exponents. The numerical estimate of g 3 is shown in Figure 7 (using the analytical results for s and s 1 s 2 ). Its decay to 0 indicates once more that s 1 s 2 s 3 does not diverge fast enough. Along the same lines, Figure 8 shows s 1 s 2 s 3 s 4 /[(x 3 0 /D) 4 T 3/2 ] together with a fit to a parabola in ln(T ). While the fit is not perfect and the data is, apart from the result from importance sampling, slightly inconsistent for large T , the suggested parabola seems plausible. The data is certainly not consistent with gap scaling.

Summary
Motivated by the question whether the sequence of the areas between trajectories of coalescing random walkers exhibits multiscaling, we have calculated the second moment and subsequently the two point correlation function exactly. The correlation function displays anti-correlations, as was discovered in the analysis of the variance of the numerical estimator of the first moment. Anti-correlations are sometimes regarded as indicating multiscaling [10,14].
Three moments are necessary to expose the absence of gap scaling, which necessitated a numerical calculation. Due to the presence of rare but important events, convergence is slow, yet an importance sampling scheme produces results consistent with logarithmic scaling, which is often observed in conjunction with multiscaling [10,14].
While the analytical proof for multiscaling in the present system is still lacking, the numerical results indicate it convincingly. Two possible avenues may be pursued in the future: An exact solution along the lines discussed in section 3.1.3, or an exact calculation based on the work by Munasinghe et al [10,14].