Test finite-ring TASEP dynamics

Can a random trajectory calculation reproduce an exact finite-state evolution? On a four-site TASEP ring, you can answer this question without a large solver: enumerate six configurations, evolve their probabilities, and compare independent trajectories at fixed times. A coordinate-Bethe mode supplies a second analytic check. Two deliberately wrong interpretations show why stationarity and a small eigenvector residual alone are insufficient tests.

Required background. Be able to construct a probability-column generator using the Markov bridge and explain the stationary current on a ring. Elementary Python helps you change the calculation; all definitions and selected outputs appear below. No quantum-mechanical preparation is needed.

Six configurations and their probability law

Use continuous-time right-moving TASEP with $N=4$ sites, $M=2$ indistinguishable particles and rate $r\gt0$ for each eligible jump. A particle moves one site to the right only if that site is empty, including the periodic bond $4\to1$. Label a configuration by its occupied sites and fix the order

$$(12,13,14,23,24,34).$$

A probability column obeys $p'(t)=Qp(t)$, with a transition rate in row destination, column source. For example, $12\to13$ is the only exit from $12$, whereas $24$ can jump to either $12$ or $34$. Thus their waiting rates are $r$ and $2r$, respectively. The diagonal entry is minus the total exit rate. This is the finite-ring convention of Golinelli and Mallick 2006, § II.A, pp. 2–3, equations (1)–(3), PDF, with an explicit rate $r$ restored.

The experiment constructs $Q$ from occupied-site bit masks and compares it with a separately specified six-state transition matrix. It checks nonnegative off-diagonal entries, zero column sums, the uniform stationary vector $\pi=\tfrac16\mathbf1$, and the exact stationary current $j=r/3$ per bond. The model fixes the regime; the convention reference distinguishes this current from the total jump rate $4r/3$.

The main time-evolution test starts from configuration $12$ with probability one, so $p(0)=e_{12}$. The input uses $r=1$ and observation times $rt=0.25,1,2$. Stationarity is a separate check, not an assumption about this initial state.

Evolve probabilities with a controlled tail

The largest exit rate is $\nu=2r$. Therefore

$$P=I+\frac{Q}{\nu}$$

is a nonnegative matrix whose columns sum to one. Expanding the exponential gives the uniformization formula

$$p(t)=e^{-\nu t}\sum_{k=0}^{\infty} \frac{(\nu t)^k}{k!}P^k p(0).$$

Each $P^kp(0)$ is a probability vector. The sum mixes discrete updates with a Poisson-distributed update count; the diagonal of $P$ supplies possible self-transitions. These self-transitions belong to the calculation, not to additional physical particle hops.

Let $w_k=e^{-\mu}\mu^k/k!$, where $\mu=\nu t$. After retaining terms through $K$, the omitted mass is $\sum_{k>K}w_k$. When $K+2\gt\mu$, successive omitted weights have ratios bounded by $\mu/(K+2)$, so

$$\sum_{k>K}w_k\leq \frac{w_{K+1}}{1-\mu/(K+2)}.$$

In exact arithmetic, this omitted mass is also the $L^1$ error of the positive truncated sum. The code stops when the upper bound is below $10^{-13}$ and does not renormalize away its missing mass. Binary64 rounding is a separate error; deterministic checks allow $2\times10^{-12}$ for it.

At $rt=0.25,1,2$, the last included degrees are $12,19,26$, with analytic tail bounds approximately $1.23\times10^{-14}$, $6.45\times10^{-14}$ and $3.54\times10^{-14}$. The method uses no numerical eigenvector decomposition as its evolution oracle.

Sample independent trajectories at fixed times

The second implementation stores a list of particle positions and constructs available moves directly. It does not read the generator matrix. With $b$ available moves, it draws a waiting time with exponential rate $rb$, chooses one of those moves uniformly, and updates the configuration. Between events the configuration remains constant.

Record the configuration at each prescribed observation time, even if no jump occurs then. The saved calculation uses 50,000 independent trajectories, NumPy’s PCG64 generator and seed 20261003, all starting at $12$. Different times on one trajectory are correlated; trajectories are independent of one another in the sampling model.

At $rt=1$, the deterministic and sampled probabilities are:

Configuration	Uniformization	Sample frequency	Pointwise 95% Wilson interval
$12$	$0.388883$	$0.387300$	$[0.383039,0.391578]$
$13$	$0.240126$	$0.236320$	$[0.232617,0.240064]$
$14$	$0.136688$	$0.139720$	$[0.136709,0.142787]$
$23$	$0.136688$	$0.137680$	$[0.134688,0.140728]$
$24$	$0.076612$	$0.077560$	$[0.075248,0.079937]$
$34$	$0.021003$	$0.021420$	$[0.020187,0.022726]$

All entries are rounded. The intervals are approximate pointwise binomial intervals. They invert the score inequality $n(\widehat p-p)^2\leq z^2p(1-p)$, with $z\approx1.96$: solving this quadratic for $p$ gives the endpoints, following Wilson 1927. The exact probabilities for $13$ and $14$ narrowly miss their intervals in this run. A nominal 95% interval does not guarantee coverage of every state at every time; these misses are retained rather than hidden by changing the seed.

For an independently justified acceptance check, if $n$ trajectories estimate one probability $p$ by $\widehat p$, apply Hoeffding 1963, p. 15, Theorem 1, equation (2.3), PDF to the Bernoulli indicators and their complements. Adding the two one-sided bounds gives

$$\Pr\!\left(|\widehat p-p|\geq\delta\right) \leq2e^{-2n\delta^2}.$$

Apply a union bound to the $m=6\times3=18$ comparisons. With $\alpha=10^{-6}$, choose

$$\delta=\sqrt{\frac{\log(2m/\alpha)}{2n}} \approx0.01319.$$

This bounds the probability of any larger sampling discrepancy by $\alpha$, under the independent-trajectory model. It requires no independence between observation times or between the six state counts. The code also allows the tiny deterministic tail and rounding error. The observed maximum absolute discrepancies are $0.0014075$, $0.0038055$ and $0.0049287$ at the three times. This conservative check is distinct from the narrower descriptive Wilson intervals. One fixed seed does not establish a Monte Carlo convergence rate.

Reconstruct a Bethe mode and a probability perturbation

The Library derivation constructs the right eigenvector

$$\psi(x,y)=z_1^xz_2^y+S z_2^xz_1^y, \qquad S=-\frac{1-z_2}{1-z_1}.$$

For $z_1=e^{2\pi i/3}$, $z_2=z_1^{-1}$, the contact and periodic conditions give an eigenvalue $E=-3r$. The vector is proportional to

$$v=(1,-2,1,1,-2,1)^{\mathsf T}.$$

The general coordinate eigenvalue and periodic equations are stated in Golinelli and Mallick 2004, § 2.2, p. 4, equations (3)–(6), PDF. The code independently reconstructs this particular mode and measures

$$\frac{\|Q\psi+3r\psi\|_2}{r\|\psi\|_2} \approx2.48\times10^{-16}.$$

It also computes the characteristic polynomial of $Q/r$ using exact rational arithmetic:

$$\det(qI-Q/r)=q(q+1)^2(q+3)(q^2+3q+4).$$

The eigenvalues are $0,-r,-r,-3r$ and $r(-3\pm i\sqrt7)/2$. The chosen mode therefore decays faster than the slowest nonstationary modes. This finite check does not prove Bethe completeness for arbitrary rings.

Because $v$ contains negative components and sums to zero, it is not a probability distribution. Instead use the separate initial condition

$$p_a(0)=\pi+a v, \qquad a=\frac1{24}.$$

It has the exact evolution $p_a(t)=\pi+a e^{-3rt}v$. Uniformization agrees with this formula to $L^1$ errors at most $6.41\times10^{-14}$ on the tested grid, consistent with the truncation bounds plus rounding. This analytic test uses a different initial state from the trajectory table.

Two mistakes that simple checks miss

Reversing the generator

On this ring $Q^{\mathsf T}\pi=0$ as well as $Q\pi=0$. The particular vector $v$ also satisfies the same left and right eigenvalue equation. Thus neither uniform stationarity nor this mode detects transposition.

The initial derivative from $12$ does:

$$\begin{aligned} Qe_{12}/r&=(-1,1,0,0,0,0)^{\mathsf T},\\ Q^{\mathsf T}e_{12}/r&=(-1,0,0,0,1,0)^{\mathsf T}. \end{aligned}$$

The first allows $12\to13$; the second puts the gain into $24$. Their $L^1$ difference is $2$. This control checks a physical transition direction, not just an abstract conservation identity.

Counting jump instants as time samples

The exit-count vector is $b=(1,2,1,1,2,1)$. With $D=\operatorname{diag}(rb)$, the embedded jump chain has transition matrix

$$J_{\mathrm{jump}}=I+QD^{-1}.$$

Its invariant distribution is proportional to the probability flux out of each state, $D\pi$, giving

$$\pi_{\mathrm{jump}}=\frac18(1,2,1,1,2,1)^{\mathsf T}.$$

This differs from the uniform time-stationary distribution by $L^1$ distance $1/3$. Faster-exiting configurations are counted more often at jump epochs. Dividing these weights by the exit rates and normalizing recovers $\pi$. These are exact stationary-distribution calculations; they have no Monte Carlo sampling uncertainty. The embedded chain is periodic, so its invariant distribution does not imply convergence from a fixed state at every successive jump number. Fixed-time trajectory sampling avoids this interpretation error.

Run and change the experiment

Download and extract the complete TASEP experiment (ZIP), then open its tasep folder. Individual files are also available: experiment.py, inputs.json, results.json, requirements.txt, and computation notes. The recorded run used Python 3.9.6, NumPy 2.0.2 and binary64 arithmetic; only NumPy is required beyond the Python standard library. With Python 3.9–3.12:

python3 -m venv .venv
.venv/bin/python -m pip install -r requirements.txt
.venv/bin/python experiment.py --check

On Windows use .venv\Scripts\python.exe. In the repository with NumPy installed, run python3 public/computations/tasep/experiment.py --check.

The check recomputes the experiment without writing files, requires identical saved inputs, compares integer counts exactly, and compares floating-point outputs with absolute tolerance $2\times10^{-12}$ and relative tolerance $10^{-9}$. Runtime versions are reported but not compared. The fixed seed supports reproduction; the separate uncertainty calculation supports statistical interpretation. Running without a flag prints fresh results; --write deliberately replaces the saved JSON. Use a separate copy for exploration. Routine site builds do not run the experiment.

Exercises

Guided: interpret a uniformization update

Starting at $12$, calculate $Pe_{12}$ for $\nu=2r$. Explain its self-transition and recover the correct derivative at $t=0$.

Hint

Use $P=I+Q/(2r)$ and expand the first two Poisson terms to first order in $t$.

Solution

$Pe_{12}=\tfrac12e_{12}+\tfrac12e_{13}$. The self-transition fills the gap between the uniformization rate $2r$ and the physical exit rate $r$. Since $e^{-2rt}=1-2rt+O(t^2)$,

$$p(t)=e_{12}+2rt(P-I)e_{12}+O(t^2) =e_{12}+rt(e_{13}-e_{12})+O(t^2).$$

Thus $p'(0)=Qe_{12}$, with no added physical jump process.

Independent: turn a signed mode into probabilities

Find every real $a$ for which $p_a(0)=\pi+av$ is a probability vector. Does $a=0.1$ work? Show that each allowed initial vector remains a probability vector under the analytic evolution.

Hint

There are only two component values: $1/6+a$ and $1/6-2a$. The components of $v$ sum to zero.

Solution

Normalization is automatic; nonnegativity requires $-1/6\leq a\leq1/12$. The choice $a=0.1$ fails because $1/6-0.2\lt0$. For $t\geq0$, multiplication by $e^{-3rt}\in(0,1]$ moves $a$ toward zero, preserving both inequalities. Hence $p_a(t)=\pi+a e^{-3rt}v$ remains nonnegative and normalized.

Transfer: remove the periodic seam

Delete the jump $4\to1$ while retaining rightward exclusion on the other bonds and no reservoirs. Starting from $12$, identify the eventual configuration. Which stationary-current, uniformity and Bethe checks must now fail?

Hint

Each successful jump increases the sum of occupied positions. This sum is bounded above, and particles cannot leave site 4.

Solution

The absorbing configuration is $34$. From every other state an allowed right jump remains, and repeated jumps eventually reach it. The stationary probability is therefore concentrated at $34$, with zero current. Uniform stationarity, the periodic formula $j=r/3$, and the old eigenmode and characteristic polynomial fail. Column sums and off-diagonal positivity still hold: the altered matrix is a valid generator for a different boundary problem. Periodic Bethe conditions cannot be carried over unchanged.

References

• Golinelli, Olivier, and Kirone Mallick. “Bethe Ansatz calculation of the spectral gap of the asymmetric exclusion process.” Journal of Physics A: Mathematical and General 37 (2004), 3321–3331. DOI: 10.1088/0305-4470/37/10/001. arXiv:cond-mat/0312371v1.
• Golinelli, Olivier, and Kirone Mallick. “The asymmetric simple exclusion process: an integrable model for non-equilibrium statistical mechanics.” Journal of Physics A: Mathematical and General 39 (2006), 12679–12705. DOI: 10.1088/0305-4470/39/41/S03. arXiv:cond-mat/0611701v1.
• Hoeffding, Wassily. “Probability Inequalities for Sums of Bounded Random Variables.” Journal of the American Statistical Association 58, no. 301 (1963), 13–30. DOI: 10.1080/01621459.1963.10500830. Open PDF.
• Wilson, Edwin B. “Probable Inference, the Law of Succession, and Statistical Inference.” Journal of the American Statistical Association 22, no. 158 (1927), 209–212. DOI: 10.1080/01621459.1927.10502953.