Introduction↩︎

The true self-avoiding walk [1] was introduced to describe a growing polymer on a lattice, where the monomer to be added to the end of the walk tries to avoid previously visited sites. It was studied using various approximate analytic tools, and it was concluded that the statistics of this walk are distinct from those of an equilibrium polymer. This self-avoiding process has found applications in several fields including chemotaxis [2], [3], but also as a tool to understand the large-scale dynamics of a family of non-reversible Monte Carlo algorithms [4]–[9]. There are also direct links to a modified “totally symmetric simple exclusion process” (TASEP) [10], [11].

Great progress in understanding this process was made in a series of papers [12]–[15] which give the full analytic solution to the motion in one spatial dimension. In particular, it was demonstrated that the motion explores a spatial extent growing in time as \(t^{2/3}\). These papers give explicit results for the distribution function of the end-to-end separation of the walk.

The founding paper [1] of this model also studied the process in the continuum with the coupled equations \[\begin{align} \frac{d {X(t)} }{d t } &=& -\frac{\partial h( {X}(t), t)}{\partial x} \tag{1} \\ \frac{dh({x},t ) }{dt} &=&\,\, \delta ( {x} - {X(t)}) \tag{2} \end{align}\] In eq. (1 ) the function \(h(x,t)\) measures the number of previous visits to the position \(X\), which acts as a repulsive potential in the motion of the end. In eq. (2 ) the repulsive potential is itself modified by visits of the walker to the position \(X(t)\). We note that we have changed notation compared to certain papers, including our own, to remain consistent with [16], the main inspiration of the present paper.

Eq. (1 ) describes the motion of a particle moving at velocity \(u=-{\partial h(t, {X}(t))}/ {\partial x}\), so it corresponds to the continuity equation \[\frac{\partial \rho}{\partial t} + \frac{\partial (u \rho)}{\partial x}=0 \label{eq:cont1}\tag{3}\] where \(\rho=\delta(x-X(t))\) is the particle density. If we now consider the spatial derivative of eq. (2 ) we find \[\frac{\partial u}{\partial t} + \frac{\partial \rho}{\partial x} =0 \label{eq:cont2}\tag{4}\] The properties of the coupled equations (3 , 4 ) for a collection of many walkers are the main subject of the present paper. Many general properties of these equations are established by [16] where these equations are shown to be the continuum limit of a discrete model, which the authors interpret as a group of bricklayers building a long wall. In order to build a uniform wall, masons move down the local gradient of height \(h\), leaving one brick on the wall each time they move. One is led to consider the collective hydrodynamic behavior of \(N\) independent bricklayers described by a continuum density distribution \(\rho(x)\), as a natural generalization of the true self-avoiding walk. The work of [16] can be placed in the context of other nonlinear models of interface growth including the model of Kardar, Parisi and Zhang (KPZ) [17] which also describes the dynamics of an interface \(h(x,t)\).

In the next section, we give the exact definition of the discrete model of bricklayers and perform simulations on this model with a variable number of bricklayers, \(N\). We consider the analytic solution of the partial differential equations for the case of bricklayers all starting together. We also study the dynamics of bricklayers starting at random positions on a periodic, circular wall.

Bricklayers↩︎

Let us recall, briefly, the discrete model of [16]. A wall is being built with bricks added on the links \(j,j+1\) of a one-dimensional lattice. The number of bricks on this link is denoted \(h_j\). Inspired by the continuum equations, we also introduce the discrete negative gradient of \(h_j\) such that \(z_j=h_{j-1} - h_j\). The bricklayers are on the nodes of the lattice with no constraints on their occupation number. At any moment, a bricklayer can jump to the right or left, leaving one brick in the appropriate column. The jumps occur at a rate \(r(z_j)\). One imposes \[\label{rcond} r(z)r(-z+1)=1,\tag{5}\] So that the local jump rate depends on the local gradient of the height of the wall, corresponding to the number of previous visits of bricklayers. In our numerical work, we choose \(r\sim\exp(\beta z)\) with \(\beta =0.4\). When a bricklayer jumps, the following changes of configuration may occur: \[(n_j,z_j),(n_{j+1},z_{j+1}) \to (n_j-1,z_j-1),(n_{j+1}+1,z_{j+1}+1)\] with rate \(n_j r(z_j)\), and \[(n_j,z_j),(n_{j-1},z_{j-1}) \to (n_j-1,z_j+1),(n_{j-1}+1,z_{j-1}-1)\] with rate \(n_j r(-z_j)\).

The number of bricklayers at site \(j\) is \(n_j\). We note that \(\sum_j n_j=N\) and \(\sum_j z_j\) are both conserved, as is the parity of \((n_j+z_j)\) for each site. In the long-time limit we expect to generate smooth distribution functions. It is natural to identify the discrete label \(j\) and the continuum variable \(x\). The discrete occupation number \(n_j\) maps to the distribution \(\rho(x)\) and \(z_j\) is linked to \(u(x)\).

We generate a starting configuration and store the randomly generated event times in a red-black tree using the C++ multiset container. This allows us to find the first event in a time that increases only logarithmically with system size. We manage the first event and calculate new event times for sites that have been modified by the configurational change. The process is repeated to evolve the system a total time \(t\).

Scaling with builder number↩︎

Here, we consider a simple generalization of previously given scaling arguments [18] for the width of the function \(\rho\) after time \(t\), starting with all builders on a single site. We consider a system with \(N\) builders, with a characteristic spatial scale which varies as \(\ell \sim N^\beta t^\alpha\) after time \(t\). The number of events is comparable to \(Nt\), which builds a wall of height \(Nt/\ell\). The gradient of height must then scale as \({ \partial h}/{\partial x} \sim N t/\ell^2\). After acting for a time \(t\) on the bricklayers this must give a scale again comparable to \(N^\beta t^\alpha = Nt^2/(N^{2\beta}t^{2\alpha})\). This implies \(\alpha=2/3\), \(\beta=1/3\), so that the extent of the distribution varies as \[\ell \sim N^{1/3} t^{2/3}\] The mean density in the occupied region then varies as \({(N/t)}^{2/3}\). We expect these results to hold when \(t\gg N\) so that the mean occupation of bricklayers in the propagating region is smaller than unity.

Simulation results↩︎

We performed simulations to study the evolution of the functions \(h(x,t)\) and \(\rho(x,t)\). We place all the bricklayers on a single starting site (the origin \(x=0\) in our figures) and at the end of the simulation record the empirical distributions. To generate high-quality data, we typically average over \(10^5\) realizations of the process. The data for different times and bricklayer number display disparate scales; when we compare the distribution functions we scale the functions to unit area and unit variance, and denote the result \(\bar \rho(\bar x,t)\) and \(\bar h(\bar x,t)\), with scaled coordinate \(\bar x\) and scaled distributions \(\bar \rho\) and \(\bar h\).

We started by studying in detail the case \(N=1\) for which explicit results are known analytically [14], confirming that the function \(\bar \rho(x,t)\) is correctly reproduced by the simulations. See the case \(N=1\) in Fig. 1. We then varied the number of bricklayers while keeping the simulation time constant at \(t=512\). We see that the double-peak structure familiar from \(N=1\) is maintained, but that the peaks are pushed to larger separations, and that the distribution of density then drops to zero more abruptly for large values of \(N\). The evolution of \(h\) is shown in Fig. 2: For small \(N\), the interface has a sharp maximum at the origin (the site where the builders were placed), but this is replaced by a broad maximum for \(N\) large.

We also studied the evolution in time when starting with \(N=512\) bricklayers, Fig. 3. At short times, \(t=2\) there is a peak at the origin, corresponding to builders that have not yet had time to move. By time \(t=32\) a dip has formed in the distribution function. The dip is shallower for the longest time shown, \(t=512\). From our simulations we were able to confirm the law in \(N^{1/3} t^{2/3}\) for the width of the distribution \(\rho(x,t)\).

Scaling solution of partial differential equations↩︎

Scaling analysis tells us that the coupled partial differential equations have a solution of the form \[\begin{align} \rho = \frac{1}{t^{2/3}} f( x/t^{2/3}) \\ u = \frac{1}{t^{1/3}} g( x/t^{2/3}) \end{align}\] These scaling relations are similar to those displayed by the solution to the KPZ equation. The partial differential equations have two conservation laws corresponding to \(\int u\,dx=0\) and \(\int \rho\, dx =N\). Substituting into the partial differential equations, we find the coupled equations of \(f\) and for \(g\) which can be written in matrix form: \[\begin{pmatrix} 3 g -2y & 3 f\\ 3 & -2y \end{pmatrix} \begin{pmatrix} f' \\ g' \end{pmatrix} = \begin{pmatrix} 2 f\\ g \end{pmatrix} \label{eq:coupled}\tag{6}\] where \(y=x/t^{2/3}\). Clearly, there is a trivial solution to these equations \(f(y)=g(y)=0\).

We used a Runge-Kutta integrator, to solve these equations from initial values of \(f\) and \(g\) at \(y=0\). We note that the value of \(f(0)\) is arbitrary, and can be chosen to be unity by rescaling of \(x\) and \(t\), and that the natural choice for a non-singular odd function \(g\) is \(g(0)=0\). In Fig. 4 we show the numerical solution of the coupled equations, eq. (6 ). We find that to high accuracy \(g(y)\) is a linear function of \(y\), while \(f(y)\) is quadratic. These results are surprising, since for large \(y\) both functions must be zero; they describe the finite time evolution of a compact distribution of bricklayers. We show below that this is possible if we introduce discontinuities in the two functions, so that the solution jumps to the trivial solution \(f(y)=g(y)=0\).

Analytic solution of the partial differential equations↩︎

The theory of hyperbolic partial differential equations is complicated by the non-uniqueness of solutions in the presence of shocks or discontinuities. The theory of the selection of the correct physical solution in such cases is subtle and will not be attempted here. We rather construct a solution from observation of the simulation data Fig. 3, informed by the numerical integration of the coupled equations Fig. 4. Empirically, we understand that the long-time limit of the function \(\rho\) is parabolic as found by the numerical integration, but it is truncated to zero at some finite \(y=y_0\), as shown by the simulations.

From the numerical integration for the scaling functions, we see that \(g= \gamma y\). If we choose \(3 \gamma =2\) the equations are simple, and we only need to solve \[3f' = 2y g' + g \label{eq:singular}\tag{7}\] for \(f\) giving $ f = A+y^2/3$. Let a discontinuity occur at \(y=y_0\), then balancing the singularities in eq. (7 ) gives \(3\Delta f-2y_0 \Delta g=0\), where \(\Delta f\) and \(\Delta g\) are the jumps in the respective functions. This allows us to relate \(A\) and \(y_0\) giving $ A= y_0^2/9\(. We also have that the total number of bricklayers is given by\)\(\int ^{y_0}_{-y_0} f(y) dy = 2 y_0 A + 2 y_0^3/9 = N\)$ so \[y_0 = {\Big ( \frac{9 N}{4} \Big )}^{1/3}\] agreeing with the previous scaling argument for the width of the distribution as a function of \(N\). We now integrate \(g(y)\) to find the wall profile: we have \(-\partial h/\partial x = (1/t^{1/3})(2/3) (x/t^{2/3}) = 2x/(3t)\) so that $ h(x,t)= B - x^2/t$ with the zero of \(h\) at \(x= t^{2/3} y_0\). Thus, \(B= (1/3) t^{1/3} y_0^2\). We find the scaling form of the solution to the partial differential equations is \[\begin{align} h(x,t)=& \frac{t^{1/3}}{3} \left ( y_0^2 - {\left (\frac{x}{t^{2/3}}\right )}^2 \right ) \Theta(y_0 -|x/t^{3/2}| )\\ \rho(x,t) =& \frac{1}{9 t^{2/3}} \left (y_0^2+ 3 {\left ( \frac{x}{t^{2/3}} \right)}^2 \right ) \Theta(y_0 -|x/t^{3/2}| )\\ u(x,t) = & \frac{2 x}{3t} \Theta(y_0 -|x/t^{3/2}| ) \label{eq:dist95rho} \end{align}\tag{8}\] The equations display discontinuous behavior at \((x/t^{2/3})=y_0\).

In Fig. 5 we plot the evolution of the empirical distribution \(\rho\) as a function of the number of bricklayers \(N\) plotting as a function of \(\bar x^2\). For these detailed comparisons, we use up to \(N=1024\) bricklayers, with simulation times up to \(t=32,384\) to ensure the convergence of the distribution functions. We see that simulations with a large number of bricklayers converge to a triangular function in this plot, indicating that the empirical distribution \(\rho(x)\) is indeed parabolic at large times. Note, however, that the numerical agreement is slightly less good for small values of \(\bar x\).

From these curves, it is also possible to measure the width of the transition region from large to small \(\rho\). An empirical fit finds this width \(w\sim t^{2/3}/N^{0.37}\). This is larger than the mean builder spacing \(t^{2/3}/N^{2/3}\). We have been unable to find an explanation for this law, which clearly goes beyond the continuum limit of the partial differential equations.

Linearized equations↩︎

When there are a large number of bricklayers distributed along the wall, and the slope of the wall is small, we can write \(\rho = \rho_0 + \delta \rho\) finding the linearized equation for the density field. \[\frac{\partial^2 \rho}{\partial t^2 } = \rho_0 \frac{\partial^2 \rho}{\partial x^2 }\] We recognize this as the wave equation with propagation speed \(c^2 = \rho_0\). The bricklayers thus move in traveling waves when they are sufficiently dense and uniform. We perform simulations by randomly placing pairs of builders (to impose even parity on \((z_i+n_i)\)) on \(L/8\) sites of a flat wall of \(L\) sites. We simulate for some time, then record the longest wavelength Fourier components of \(\rho\). In Fig. 6 we plot the autocorrelation of the resulting time series. We see indeed coherent oscillations, corresponding to two counter-rotating waves on the ring.

Conclusions↩︎

We have found a (possibly non-unique) analytic solution to the partial differential equation describing a generalized multi-walker true self-avoiding motion. This solution is compatible with the numerical solution we find to the model of [16]. Both the interface profile \(h\) and the bricklayer density \(\rho\) are simple quadratic functions of the position, but are truncated at some finite distance.

The true self-avoiding walk has been shown to have application in describing the dynamics of non-reversible Monte Carlo simulations. A multi-agent (parallel) version of this algorithm has been described in [19]. It will be interesting to see if the ideas of the present paper can also be applied to such non-reversible Monte Carlo algorithms in two or more dimensions with multiple agents. We note also that with certain parameter values non-reversible Monte Carlo methods can generate oscillating states, where the convergence of thermodynamic averages becomes difficult to control. The oscillating states found in our simulations of the bricklayer model are perhaps linked to this point.

References↩︎