Exact solution of two interacting run-and-tumble random walkers with finite tumble duration

We now precisely define our lattice model of two run-and-tumble random walkers. We consider two particles moving under stochastic dynamics on a periodic one-dimensional lattice of L sites. Each particle occupies one lattice site and has an internal velocity state $\sigma_i= 0, +1, -1$. A value $\sigma = \pm 1$ (hereafter denoted simply $+$ or  −) indicates a direction of motion to the right or left respectively; a value $\newcommand{\tum}{{0}} \sigma_{i}=\tum$ indicates that the particle is in a tumbling state and remains stationary on its site. Due to the translational invariance of the system, a microscopic configuration is fully specified by $1\leqslant n < L$, the distance between the two particles in units of the lattice spacing, and the two particle velocities, $\sigma_1$ and $\sigma_2$. A right-moving particle ($\sigma_i=+$) hops one site to the right with rate γ ; likewise, a left-moving particle ($\sigma_i=-$) hops with rate γ to the left. However when the target site is occupied by another particle, hopping is not allowed due to the hard-core exclusion interaction. Such random walks with non-intersecting paths, which implement the hard-core exclusion interaction, have a long history in physics [31].

When a particle enters a tumbling state $\sigma_i=0$, the particle stops hopping. The run lengths and tumble durations are both Poisson-distributed. The particle enters the tumbling state from a running state with rate $\tilde \alpha$ and re-enters a running state from the tumbling state with rate $\tilde \beta$. In the following we shall consider tumbling rates scaled by the hopping rate γ: $\alpha = \tilde \alpha/\gamma$ and $\beta = \tilde \beta/\gamma$. A spatiotemporal plot of a simulation of the model is provided in figure 1.

Zoom In Zoom Out Reset image size

Our main result is to obtain the exact steady state, which is unique due to the finite system size. As already mentioned, each configuration of our model is uniquely described by the distance, n, between the two particles and their respective internal velocity states $\sigma_{1}, \sigma_{2}$. The exact form for the probability distribution in the steady state, $P_{\sigma_{1}\sigma_{2}}(n)$, is

$$\begin{align} \displaystyle P_{\sigma_{1}\sigma_{2}}(n) &= a_{\sigma_{1}\sigma_{2}}(1) + a_{\sigma_{1}\sigma_{2}}(z_+) z_{+}^{-n+1} + a_{\sigma_{1}\sigma_{2}}(1/z_+) z_{+}^{n-1} \nonumber \\ &\quad + a_{\sigma_{1}\sigma_{2}}(z_-) z_{-}^{-n+1} + a_{\sigma_{1}\sigma_{2}}(1/z_-) z_{-}^{n-1} + w_{\sigma_{1}\sigma_{2}}^{(0)} \delta_{n,1} + w_{\sigma_{1}\sigma_{2}}^{(1)} \delta_{n,L-1}, \label{P intro} \nonumber \end{align} \tag{ 1 }$$

where the amplitudes $a_{\sigma_{1}\sigma_{2}}(z)$, $w_{\sigma_{1}\sigma_{2}}^{(0)}$ and $w_{\sigma_{1}\sigma_{2}}^{(1)}$, the factors $z_{+}$ and $z_{-}$ are functions of the model parameters α, β ($z_{+}$ and $z_{-}$ are implicitly defined in (35) and (36), but we leave the full expressions to the supplementary material (stacks.iop.org/JPhysA/50/375601/mmedia) in lines 191 and 192) and L and $\delta_{n, m}$ is the Kronecker delta symbol. In other words, this distribution comprises a constant part and terms that vary exponentially with the particle separation n with further contributions in states where the two particles are next to each other ($n=1$ or $n=L-1$).

We may understand this distribution by considering the dynamics of the jamming that occurs between the interacting run-and-tumble particles. Two particles with equal and opposite velocities collide so that the two particles are on neighbouring sites in either the $(+-, n=1)$ configuration or its symmetrically related counterpart $(-+, n=L-1)$. This results in a microscopically jammed configuration as neither particle can hop freely until one of them changes orientation. Furthermore, the system cannot change its configuration until one particle starts tumbling. This waiting time is reflected by a delta symbol in the probability of these jammed configurations i.e. $w_{\sigma_{1}\sigma_{2}}^{(0)}$ is nonzero in $P_{+-}(1)$, and $w_{-+}^{(1)}$ is nonzero in $P_{-+}(L-1)$. The second part of the interaction involves unjamming. Eventually the system enters a jammed configuration of type $\newcommand{\tum}{{0}} (+\tum, n=1)$ or one of the symmetric counterparts, in which one of the two particles involved in the collision has begun tumbling. These configurations therefore also each contain a delta symbol in their probabilities. There is also an enhanced probability of entering the $\newcommand{\tum}{{0}} (\tum \tum, n=1)$ state in which both adjacent particles are tumbling from these jammed tumbling configurations, which in turn generates delta symbols in $\newcommand{\tum}{{0}} P_{\tum \tum}(1)$, so that $\newcommand{\tum}{{0}} w_{\tum\tum}^{(0)}$ and $\newcommand{\tum}{{0}} w_{\tum\tum}^{(1)}$ are non zero. On the other hand $w_{++}^{(0, 1)}$ and $w_{--}^{(0, 1)}$ are each zero, as there are no delta-symbol contributions to the probabilities in these velocity sectors.

Particles unjam by leaving a jammed tumbling configuration; that is when the tumbling particle exits tumbling with an orientation different from the one it had on collision. At this point both particles move in the same direction. Due to the stochasticity of the hopping between lattice sites, some broadening of the separation between the two particles will occur, even though they started off next to each other. This broadening generates a spatially decaying component in the probability with exponential decay length scale $\xi_{+} = 1/\ln z_{+}(\alpha, \beta)$ (analogous to ξ in [11]), that is a function of the tumbling rates α and β, and is apparent in all the velocity sectors. Note that in an equilibrium system, an exponentially decaying probability arises from a linear potential with a positive gradient, or—equivalently—a constant attractive force. In this nonequilibrium system, such forces emerge from irreversible collisions between particles.

Finally, if one of the particles enters a tumbling state when the particles are separated and freely moving, the result is a configuration where one particle is stationary and tumbling and the other is hopping freely. The freely moving particle either hops towards or away from the tumbling particle. This contribution to the stationary probability distribution is characterised by an exponential decay, but with a new length scale $\xi_{-} = 1/ \ln z_{-}(\alpha, \beta)$ for which there is no equivalent in [11]. This completes our discussion of the different microscopic mechanisms that lead to (1).

As noted in the introduction, this lattice-based model is an approximation to the real-world situation of continuous space and time. In order to recover continuum dynamics, we take the lattice spacing as $\ell/L$ where $\ell$ is the physical system size and let $L \to \infty$. In order to keep the physical velocity

$$\begin{align} \displaystyle v = \gamma \frac{\ell}{L} \label{vdef} \nonumber \end{align} \tag{ 2 }$$

invariant we also scale the hopping rate γ with system size

$$\begin{align} \displaystyle \gamma = L/\ell. \label{gscale} \nonumber \end{align} \tag{ 3 }$$

Then $v = 1$ and the scaled tumbling rates (ratio of bare tumbling rates to hopping rates) scale as $1/L$:

$$\begin{align} \displaystyle \alpha = \frac{\tilde \alpha}{\gamma} = \frac{\phi}{L} \quad \beta = \frac{\tilde \beta}{\gamma} = \frac{\theta}{L} \label{scaling limit def} \nonumber \end{align} \tag{ 4 }$$

where $\phi = \ell \tilde \alpha$ and $\theta= \ell \tilde \beta$ are dimensionless constants. Thus in this scaling limit the particles undergo a ballistic motion with velocity $v=1$ interrupted by collision events and stochastic tumbles.

In the scaling limit, the steady-state probabilities of walkers at a separation y have the following form:

$$\begin{align} \displaystyle P_{++}(y) = a_{++} + b_{++}[\delta(y) + \delta(\ell-y)] + c_{++}[{\rm e}^{- y/\kappa} + {\rm e}^{-(\ell-y)/\kappa}] \label{++ scaling limit} \nonumber \end{align} \tag{ 5 }$$

$$\begin{align} \displaystyle P_{+-}(y) = a_{+-} + c^{(0)}_{+-} {\rm e}^{- y/\kappa} + c^{(1)}_{+-} {\rm e}^{- (\ell-y)/\kappa} + w_{+-}\delta(y) \nonumber \end{align} \tag{ 6 }$$

$$\begin{align} \newcommand{\tum}{{0}} \displaystyle P_{+\tum}(y) = a_{+\tum} + c^{(0)}_{+\tum} {\rm e}^{- y/\kappa} + c^{(1)}_{+\tum}{\rm e}^{- (\ell-y)/\kappa} + w_{+\tum}\delta(y) \nonumber \end{align} \tag{ 7 }$$

$$\begin{align} \newcommand{\tum}{{0}} \displaystyle P_{\tum \tum}(y) = a_{\tum \tum} + c_{\tum \tum}[{\rm e}^{- y/\kappa} + {\rm e}^{-(\ell - y)/\kappa}] + w_{\tum \tum}[\delta(y) + \delta(\ell - y)] \label{tum tum scaling limit} \nonumber \end{align} \tag{ 8 }$$

where the length scale κ is given by

$$\begin{align} \displaystyle \label{xi} \kappa = \ell \left(\frac{2}{(\theta +\phi) (\theta +2 \phi)}\right)^{1/2}. \nonumber \end{align} \tag{ 9 }$$

The amplitudes $a_{\sigma_{1}\sigma_{2}}$, $b_{\sigma_{1}\sigma_{2}}$ and $c_{\sigma_{1}\sigma_{2}}$ derive from the amplitudes $a_{\sigma_{1}\sigma_{2}}(z)$ that appear in (1). They are functions of the dimensionless parameters θ and ϕ, and are specified explicitly in section 4. Superscripts appear where these amplitudes are different at separation $y=0$ and $y=\ell$. The scaling limit found in [11] may be recovered by taking $\tilde{\beta} \rightarrow \infty$ for (5)–(8), as we show in section 4.5.

A feature of this scaling limit is that in the $++$ (and symmetric −−) sector the terms in (1) containing $z_{+}$ to some power have become delta functions. Therefore they have gone from a finite length scale $\xi_{+} = 1/ \ln z_{+}$ to a delta-function one. The origin of this vanishing lengthscale lies in the fact that the runs are no longer described by stochastic hops, which in the lattice model led to broadening of the particle separation. In the other sectors the $z_{+}$ terms disappear for the same reason. We also saw this for the analogous lengthscale ξ in [11]. The second length scale $\xi_{-} = 1/\ln z_{-}$, however, remains finite and present in all velocity sectors in the scaling limit resulting in the decay length $\kappa = \frac{\ell}{L} \frac{1}{\ln z_-}$ (9). This is because the tumble duration—and hence, distance travelled by a moving particle when the other particle is tumbling—remains finite in this limit.

Equations (5)–(8) demonstrate the rich structure that nonequilibrium steady states may have in comparison to their equilibrium counterparts. The rest of this paper sets out the derivation of these results.

To progress towards (31), we now note that the determinant of A, (34), can be written as the following polynomial fraction

$$\begin{align} \displaystyle \frac{\det A(x)}{1-x} = \frac{k}{x^{3}}q(x) \label{det A factorised} \nonumber \end{align} \tag{ 35 }$$

where

$$\begin{align} \displaystyle q(x) \equiv (x-1)(x-z_{+})(x-1/z_{+})(x-z_{-})(x-1/z_{-}), \label{qdef} \nonumber \end{align} \tag{ 36 }$$

and

$$\begin{align} \displaystyle k =\frac{\beta}{2} (2+\alpha)(1+\alpha +\beta)\label{kdef} \nonumber \end{align} \tag{ 37 }$$

is a constant. In expression (36), $z_{+}$ and $z_{-}$ are independent roots of the determinant, and $1/z_{+}$ and $1/z_{-}$ are their inverses. There is also a root at $z=1$. This furnishes the five roots $z_\rho$ that appear in (31). To be explicit: $z_0 = 1$, $z_1 = z_{+}$, $z_2=1/z_{+}$, $z_3=z_{-}$ and $z_4=1/z_{-}$.

We now manipulate the expression in (33) with a view to writing it in the form of (31). Given that we may rewrite $\det A$ as a polynomial fraction (35, 36), we write (33) as

$$\begin{align} \newcommand{\adj}{{{\rm adj}}} \displaystyle G_{\sigma_{1}\sigma_{2}} &= \frac{-x^{3}(1-x)\sum\nolimits_{j} \adj A_{\sigma_{1}\sigma_{2},j}(x)b_{j}(x)}{k(x-1)^2(x-z_{+})(x-1/z_{+})(x-z_{-})(x-1/z_{-})} \nonumber \\ &= \frac{x^{3}\sum\nolimits_{j} \adj A_{\sigma_{1}\sigma_{2},j}(x)b_{j}(x)}{kq(x)} \cdot \label{Gadic} \nonumber \end{align} \tag{ 38 }$$

We may separate terms in $G_{\sigma_{1}\sigma_{2}}(x)$ as follows

$$\begin{align} \newcommand{\adj}{{{\rm adj}}} \displaystyle G_{\sigma_{1}\sigma_{2}} = \frac{x^{3}\sum\nolimits_{j} \adj A_{\sigma_{1}\sigma_{2},j}(x)b_{j}(x)}{kq(x)} = \frac{x p_{\sigma_{1}\sigma_{2}}(x)}{q(x)} + \frac{\tilde H_{\sigma_{1}\sigma_{2}}(x)}{q(x)}, \label{GF general} \nonumber \end{align} \tag{ 39 }$$

where each combination $x p_{\sigma_{1}\sigma_{2}}(x)$ is a polynomial of degree less than $x^{L}$ and $\tilde H_{\sigma_{1}\sigma_{2}}(x)$ is a polynomial with a lowest order term $x^{L}$. Thus $x p_{\sigma_{1}\sigma_{2}}(x)$ contains all the terms of $\newcommand{\adj}{{{\rm adj}}} x^{3}\sum\nolimits_{j} \adj A_{\sigma_{1}\sigma_{2}, j}(x)b_{j}(x)$ with degree less than L.

We now show that the only terms in $xp_{\sigma_{1}\sigma_{2}}$ of order greater than $x^{5}$ are of order $x^{6}$ and $x^{L-1}$. To do this, we consider those terms in $\newcommand{\adj}{{{\rm adj}}} x^{3}\adj A_{\sigma_{1}\sigma_{2}, j}(x)b_{j}(x)$ that are $O (x^{m})$ where $5 < m < L-1$. Cramer's rule allows us to write

$$\begin{align} \newcommand{\adj}{{{\rm adj}}} \displaystyle \adj A_{i,j}(x)b_{j}(x)= \det A_{i}, \nonumber \end{align} \tag{ 40 }$$

where $A_{i}$ is the matrix formed by replacing the ith column of A with $\underline{b}$. Then we see that $O(x^{6})$ terms in $x^3\det A_{i}$ must come from $\mu^{3}$ terms in $\det A_i$. Likewise, the $O(x^{L-1})$ terms in $x^3\det A_{i}$ must come from multiplying $x^{L-1}$ terms in $\underline{b}$ by $\nu(x){\hspace{0pt}}^{3}$ terms in $\det A_{i}$. Since $\underline{b}$ only contains terms of order $O(1)$ and $O(x^{L-1})$ one can check that all other terms in $x^3\det A_{i}$ are $O(x^m)$ where either $m > L-1$ or $m< 6$.

We separate each $p_{\sigma_{1}\sigma_{2}}(x)$ into those terms that will allow partial fraction decomposition, and those that do not:

$$\begin{align} \displaystyle \frac{xp_{\sigma_{1}\sigma_{2}}(x)}{q(x)} + \frac{\tilde H_{\sigma_{1}\sigma_{2}}(x)}{q(x)} = x\frac{J_{\sigma_{1}\sigma_{2}}(x)}{q(x)} + \frac{K_{\sigma_{1}\sigma_{2}}(x)}{q(x)} + \frac{\tilde H_{\sigma_{1}\sigma_{2}}^{\prime}(x)}{q(x)}, \label{powersep} \nonumber \end{align} \tag{ 41 }$$

where $J_{\sigma_{1}\sigma_{2}}(x)$ takes terms from $p_{\sigma_{1}\sigma_{2}}(x)$ amenable to partial fraction decomposition (ie. those of order less than $q(x)$) and is therefore a polynomial of order $x^{4}$ or less, and $K_{\sigma_{1}\sigma_{2}}(x)$ takes the higher order terms from $p_{\sigma_{1}\sigma_{2}}(x)$. However, we desire an expression $K_{\sigma_{1}\sigma_{2}}(x)/q(x)$ that can be cast as a polynomial rather than a rational function so that we can read off its contribution to the probability. To this end, we define

$$\begin{align} \displaystyle K_{\sigma_{1}\sigma_{2}}(x)/q(x) \equiv w_{\sigma_{1}\sigma_{2}}^{(0)}x + w_{\sigma_{1}\sigma_{2}}^{(1)}x^{L-1}, \nonumber \end{align} \tag{ 42 }$$

where $w_{\sigma_{1}\sigma_{2}}^{(0)}$ is equal to the ratio of the coefficient of the $x^{6}$ term in $xp_{\sigma_{1}\sigma_{2}}(x)$ with the coefficient of the $x^{0}$ term in $q(x)$ and $w_{\sigma_{1}\sigma_{2}}^{(1)}$ is equal to the ratio of the coefficient of the $x^{L-1}$ term in $xp_{\sigma_{1}\sigma_{2}}(x)$ with the coefficient of the $x^{5}$ term in $q(x)$. In order to factorise $K_{\sigma_{1}\sigma_{2}}(x)$ by $q(x)$, we add to $K_{\sigma_{1}\sigma_{2}}(x)$ any terms required, in addition to the $x^{6}$ and $x^{L-1}$ terms in $xp_{\sigma_{1}\sigma_{2}}(x)$ already present. If these added terms are of degree less than 5, then we subtract them from $p_{\sigma_{1}\sigma_{2}}(x)$. On the other hand, if the added terms are of degree greater than $L-1$ (recall, we have already shown that there are no further terms between $x^{5}$ and $x^{L}$), we subtract them from $\tilde H_{\sigma_{1}\sigma_{2}}(x)$, which in turn becomes $\tilde H_{\sigma_{1}\sigma_{2}}^{\prime}(x)$.

We now return to our expressions $J_{\sigma_{1}\sigma_{2}}(x)$ in (41), which are amenable to partial fraction decomposition. A remarkable simplification occurs when we use Heaviside's 'cover-up' method for the partial-fraction expansion of a rational function [32], on the fraction $\frac{J_{\sigma_{1}\sigma_{2}}(x)}{q(x)}$. The method may be used whenever the denominator of a rational fraction can be factorised into distinct linear factors. As $q(x)$ can be written in this form, and that each $J_{\sigma_{1}\sigma_{2}}(x)$ is a polynomial, and therefore the method can be applied to our fraction, which yields

$$\begin{align} \displaystyle \frac{J_{\sigma_{1}\sigma_{2}}(x)}{q(x)} &= \frac{J_{\sigma_{1}\sigma_{2}}(x)}{(x-z_{1})(x-z_{2})...(x-z_{n})} \nonumber \\ &= \frac{J_{\sigma_{1}\sigma_{2}}(z_{1})}{(z_{1}-z_{2})...(z_{1}-z_{n})} \cdot \frac{1}{x-z_{1}} +... + \frac{J_{\sigma_{1}\sigma_{2}}(z_{n})}{(z_{n}-z_{1})...(z_{n}-z_{n-1})} \cdot \frac{1}{x-z_{n}} \nonumber \end{align} \tag{ 43 }$$

where $z_{1},..., z_{n}$ are the roots of $q(x)$. The numerators, $J_{\sigma_{1}\sigma_{2}}(z_{i})$, are found by covering up the factor $x-z_{i}$ in $\frac{J_{\sigma_{1}\sigma_{2}}(x)}{q(x)}$, and setting $x=z_{i}$ in the rest of the expression. The terms involving $J_{\sigma_{1}\sigma_{2}}$ are now in the form of $\frac{xM_{\sigma_{1}\sigma_{2}}(z_{\rho})}{(1-z_{\rho}x)}$ of (31) and so straightforwardly invertible. We can ignore the expressions within $\tilde H_{\sigma_{1}\sigma_{2}}^{\prime}(x)$ entirely as they do not contribute. We may therefore write the generating function in general as

$$\begin{align} \displaystyle G_{\sigma_{1}\sigma_{2}} = \sum_{\rho} \left[\frac{xJ_{\sigma_{1}\sigma_{2}}(z_{\rho})}{[q(x)/(x-z_{\rho})]_{|x=z_{\rho}}}\frac{1}{(x-z_{\rho})} \right] + w_{\sigma_{1}\sigma_{2}}^{(0)}x + w_{\sigma_{1}\sigma_{2}}^{(1)}x^{L-1} + x^{L}\tilde H^{\prime}_{\sigma_{1}\sigma_{2}}. \label{general GF form} \nonumber \end{align} \tag{ 44 }$$

We then write an expression for the steady-state probabilities of the form in (32)

$$\begin{align} \displaystyle P_{\sigma_{1}\sigma_{2}}(n) = \sum_{\rho}a_{\sigma_{1}\sigma_{2}}(z_{\rho}) z_{\rho}^{-n+1} + w_{\sigma_{1}\sigma_{2}}^{(0)} \delta_{n,1} + w_{\sigma_{1}\sigma_{2}}^{(1)} \delta_{L-1,1}, \label{specific P} \nonumber \end{align} \tag{ 45 }$$

where

$$\begin{align} \displaystyle a_{\sigma_{1}\sigma_{2}}(z_{\rho}) = \frac{-J_{\sigma_{1}\sigma_{2}}(z_{\rho})}{[q(x)/(x-z_{\rho})]_{|x=z_{j=\rho}}} \cdot \label{define a} \nonumber \end{align} \tag{ 46 }$$

It remains to determine which weights $w_{\sigma_{1}\sigma_{2}}^{(i)}$ are non-zero in their corresponding velocity sectors $\sigma_{1}\sigma_{2}$. We proceed column-by-column in A, replacing each with $\underline{b}$. Thanks to the symmetries $G_{+-}(x) = G_{-+}(L-x)$ and $\newcommand{\tum}{{0}} G_{+\tum}(x) = G_{\tum+}(L-x)$, we are only required to solve for four generating functions $\newcommand{\tum}{{0}} G_{++}(x), G_{+-}(x), G_{+\tum}(x)$ and $\newcommand{\tum}{{0}} G_{\tum \tum}(x)$. For $G_{++}$, as a μ is eliminated (replaced by $b_{1}$), it is not possible to get a $O(x^{6})$ term, nor a $O(x^{L-1})$ term as a ν is also eliminated. For $G_{+-}(x)$, we have sufficient factors of μ to get an $O(x^{6})$ term but no diagonal $x^{L-1}$ terms for $O(x^{L-1})$ terms. For $\newcommand{\tum}{{0}} G_{+ \tum}(x)$, we get an $O(x^{L-1})$ term only for the similar reasons. For $\newcommand{\tum}{{0}} G_{\tum \tum}$ both $\mu^{3}$ and $\nu^{3}$ terms are possible, and so $\newcommand{\tum}{{0}} G_{\tum \tum}(x)$ can possess both $O(x^{6})$ and $O(x^{L-1})$ terms.

We complete our derivation by briefly describing how to determine the constants $P_{++}(1), P_{+-}(1)$ and $\newcommand{\tum}{{0}} P_{+\tum}(1)$. We find two of these as yet undetermined constants by imposing the condition that the generating functions must not diverge at any x. As $G_{i}$ has poles at each of the roots, this condition implies that the numerator, $\newcommand{\adj}{{{\rm adj}}} \adj A_{i, j}(x)\tilde{b}_{j}(x)$, has to cancel the determinant poles and thus must equal zero at all of the roots. At $x=1$, the numerator is automatically zero. It remains to impose pole cancellation for the roots $z=z_{+}, 1/z_{+}, z_{-}, 1/z_{-}$. Although there are 24 simultaneous equations following from this condition, we find that only two are linearly independent. Therefore we can find any two of $P_{++}(1), P_{+-}(1)$ and $\newcommand{\tum}{{0}} P_{+\tum}(1)$. The remaining constant is found by imposing normalisation: $\sum\nolimits_{\sigma_{1}\sigma_{2}}\sum\nolimits_{n=1}^{L-1}P_{\sigma_{1}\sigma_{2}}(n)=1$.

We have derived the general form of the steady-state probability distribution, (45) and (46). However, there remain a number of expressions that we have not presented explicitly in terms of the model parameters α and β, namely the roots $z_{\rho} = z_{\rho}[\alpha, \beta]$, the weights $w_{\sigma_{1}\sigma_{2}}=w_{\sigma_{1}\sigma_{2}}[\alpha, \beta]$, the amplitudes $a_{\sigma_{1}\sigma_{2}}(z_{\rho})=a_{\sigma_{1}\sigma_{2}}(z_{\rho})[\alpha, \beta]$, and the constants $P_{++}(1)=P_{++}(1)[\alpha, \beta]$, $P_{+-}(1)=P_{+-}(1)[\alpha, \beta]$ and $\newcommand{\tum}{{0}} P_{+\tum}(1)=P_{+\tum}(1)[\alpha, \beta]$. It is possible to find these expressions explicitly, but due to their unwieldy form we consign the details to a Mathematica notebook in the supplementary material. The notebook performs an exact analytic calculation of the probability distribution up to the normalisation of the distributions $P_{\sigma_{1}\sigma_{2}}(n)$. Normalisation for a specific set of model parameters is achieved numerically.

A comparison of a simulation with our analytic solution for particular values of the model parameters is shown in figure 2, showing complete agreement. As in [11], we present the results in the form of effective potentials, $V(x) = -\ln P(x)$. For simplicity, we plot only the four independent sectors, in which the particles are approaching ($+-$ and $\newcommand{\tum}{{0}} +\tum$ sectors) or maintain a constant (average) separation ($++$ and $\newcommand{\tum}{{0}} \tum\tum$ sectors). We see there is an attraction towards low separations, $n\ll L$, in the sectors where the particles are approaching, and that the characteristic lengthscales differ between these two approaching sectors.

Zoom In Zoom Out Reset image size

Since each generating function $G_{\sigma_{1}\sigma_{2}}(x)$ is a finite sum (see (16)), it cannot diverge at any x. Consequently each pole in (38), corresponding to a root $z_{\rho}$ of $\det A$, must be cancelled by a zero in the numerator. Each such condition leads to a linear equation in $P_{++}(1)$, $P_{+-}(1)$ and $\newcommand{\tum}{{0}} P_{+\tum}(1)$. As previously noted, there are only two linearly independent equations in these quantities, with a further condition (namely, normalisation) required to determine them all.

The numerator of (38) is always zero at the root $x=1$, which does not provide any information about $P_{++}(1), P_{+-}(1)$ and $\newcommand{\tum}{{0}} P_{+\tum}(1)$. At the other roots $x=z_{+}, 1/z_{+}, z_{-}, 1/z_{-}$, we impose the condition

$$\begin{align} \newcommand{\adj}{{{\rm adj}}} \displaystyle \sum_{j} \adj A_{\sigma_{1}\sigma_{2},j}(z_{\rho})\,b_{j}(z_{\rho}) = 0. \label{pole cancellation condition} \nonumber \end{align} \tag{ 48 }$$

We find the two desired linearly independent conditions by taking $z=z_{+}$ and $z=z_{-}$ in (48) with $\sigma_{1}\sigma_{2}=++$. Using the expression (28) for $b_{j}(x) = 0$, each of these conditions takes the form

$$\begin{align} \newcommand{\tum}{{0}} \displaystyle &\tilde{A}_{++,1}(z_{\rho})P_{++}(1) + \tilde{A}_{++,2}(z_{\rho})P_{+-}(1) + \tilde{A}_{++,5}(z_{\rho})P_{+ \tum}(1) \nonumber \\ &\qquad \qquad = z_{\rho}^{L-1} \left[\tilde{A}_{++,1}(z_{\rho})P_{++}(1) + \tilde{A}_{++,3}(z_{\rho})P_{+-}(1) + \tilde{A}_{++,4}(z_{\rho})P_{+ \tum}(1) \right] \label{root cancelling G++ at z} \nonumber \end{align} \tag{ 49 }$$

where, for convenience, we introduce the notation $\newcommand{\adj}{{{\rm adj}}} \tilde{A}_{\sigma_{1}\sigma_{2}, j}(x) \equiv \adj A_{\sigma_{1}\sigma_{2}, j}$.

To apply these two conditions, we need to know the location of the roots $z_{+}$ and $z_{-}$ of $\det A$, as defined by (35). By expanding $z_{\pm}$ about 1 in powers of $1/\sqrt{L}$ in the explicit expression (34) for the determinant, one finds that

$$\begin{align} \displaystyle \label{z+sl} z_{+} \sim 1 + \frac{\sqrt{2\phi}}{\sqrt{L}} +\frac{\phi}{L} +O\left(L^{-3/2}\right) \nonumber \end{align} \tag{ 50 }$$

$$\begin{align} \displaystyle z_{-} \sim 1 + \sqrt{\frac{(\theta +\phi) (\theta +2 \phi)}{2}} \frac{1}{L} +O\left(L^{-3/2}\right). \nonumber \end{align} \tag{ 51 }$$

When we substitute these roots into (49), we find that $z_{+}^{L-1}\to\infty$, and so the terms in square brackets on the right-hand side of (49) need to cancel at this root. At $z_{-}$, the factor $z_{-}^{L-1}$ approaches ${\rm e}^\lambda$ where $\lambda = \lim\nolimits_{L \to\infty} L [z_{-} - 1]$ is given by

$$\begin{align} \displaystyle \lambda = \left(\frac{(\theta +\phi) (\theta +2 \phi)}{2}\right)^{1/2}. \nonumber \end{align} \tag{ 52 }$$

For future reference, it is also helpful to note the locations of the reciprocal roots

$$\begin{align} \displaystyle 1/z_{+} \sim 1 - \frac{\sqrt{2\phi}}{\sqrt{L}} +\frac{\phi}{L} +O\left(L^{-3/2}\right) \nonumber \end{align} \tag{ 53 }$$

$$\begin{align} \displaystyle 1/z_{-} \sim 1 - \sqrt{\frac{(\theta +\phi) (\theta +2 \phi)}{2}} \frac{1}{L} +O\left(L^{-3/2}\right). \nonumber \end{align} \tag{ 54 }$$

The next step is to determine the leading large-L forms of the adjugate elements $\tilde{A}_{\sigma_{1}\sigma_{2}, j}(x)$ appearing in (49) at each of the roots. All subleading terms will vanish in the scaling limit. To identify these leading terms, we require explicit expressions for $\tilde{A}_{\sigma_{1}\sigma_{2}, j}(x)$ in the constants α, β and the functions $\mu(x) = x-(1+\alpha)$ and $\nu(x) = x^{-1}-(1+\alpha) = \mu(x^{-1})$. These can be obtained most straightforwardly using a computational algebra package such as Mathematica. For example, one finds

$$\begin{align} \displaystyle \tilde{A}_{++,1}(x) &= -\frac{1}{4} \beta \bigg(\beta ^2 (\alpha +2 \mu) (\alpha +2 \nu)-\beta (\alpha +2 \mu) (\mu +\nu) (\alpha +2 \nu) \nonumber \\ &\quad {}+ 2 \mu \nu (\alpha (\mu +\nu)+2 \mu \nu)\bigg) \label{A++,1(zpl)} \nonumber \end{align} \tag{ 55 }$$

where μ and ν, defined in (29), are functions of x. Substituting the L-dependent expressions for α and β, (47), along with the large L form of $z_+$ (50) into the above expression yields the large-L result

$$\begin{align} \displaystyle \tilde{A}_{++,1}(z_{+}) \sim -\frac{4\theta \phi^{2}}{L^{3}} + O\left(L^{-7/2} \right). \label{example A} \nonumber \end{align} \tag{ 56 }$$

Using the same method, one can find the leading terms of each of the adjugate elements in (49) at each of the roots $z_{\rho}$. The results are summarised in table 1.

Table 1. Adjugate elements in the scaling limit required to evaluate $P_{++}(1)$ and $P_{+-}(1)$ and $J_{++}$.

	$x=1$	$x=z_{+}$	$x=1/z_{-}$
$\tilde{A}_{++, 1}$	$-\frac{\theta ^2 \phi ^2 \zeta}{4 L^5}$	$-\frac{4 \theta \phi ^2}{L^3}$	$\frac{\theta ^3 \zeta^2}{4 L^5}$
$\tilde{A}_{++, 2}$	$-\frac{\theta ^2 \phi ^2 \zeta}{4 L^5}$	$-\frac{\theta ^2 \phi ^2}{L^4}$	$-\frac{\theta ^2 \phi \zeta \left(2\lambda+\zeta \right)}{4 L^5}$
$\tilde{A}_{++, 3}$	$-\frac{\theta ^2 \phi ^2 \zeta}{4 L^5}$	$-\frac{\theta ^2 \phi ^2}{L^4}$	$-\frac{\theta ^2 \phi \zeta \left(-2\lambda+\zeta \right)}{4 L^5}$
$\tilde{A}_{++, 4}$	$-\frac{\theta ^2 \phi ^2 \zeta}{4 L^5}$	$\frac{\sqrt{2} \theta ^2 \phi ^{3/2}}{L^{7/2}}$	$-\frac{\theta ^3 \zeta \left(2\lambda-2\zeta \right)}{8 L^5}$
$\tilde{A}_{++, 5}$	$-\frac{\theta ^2 \phi ^2 \zeta}{4 L^5}$	$-\frac{\sqrt{2} \theta ^2 \phi ^{3/2}}{L^{7/2}}$	$\frac{\theta ^3 \zeta \left(2\lambda+2 \zeta \right)}{8 L^5}$

Now, solving the two equations arising from substituting $x=z_+$ and $x=1/z_-$ into (49) we find for large L

$$\begin{align} \newcommand{\tum}{{0}} \displaystyle P_{++}(1) \sim \frac{\theta}{2 \sqrt{2L\phi}} P_{+\tum}(1) \label{P++(1)} \nonumber \end{align} \tag{ 57 }$$

$$\begin{align} \newcommand{\tum}{{0}} \displaystyle P_{+-}(1) \sim \frac{\theta \left({\rm e}^{\lambda} (\zeta +\lambda)-\zeta +\lambda \right)}{\phi \left(\zeta \left({\rm e}^{\lambda}-1\right)+2 \left({\rm e}^{\lambda}+1\right) \lambda \right)} P_{+\tum}(1) \label{P+-(1)} \nonumber \end{align} \tag{ 58 }$$

where

$$\begin{align} \displaystyle \zeta \equiv \theta + 2\phi \label{define zeta} \nonumber \end{align} \tag{ 59 }$$

$$\begin{align} \displaystyle \eta \equiv \theta + \phi. \label{define eta} \nonumber \end{align} \tag{ 60 }$$

The remaining constant $\newcommand{\tum}{{0}} P_{+\tum}(1)$ will be found by normalisation (see section 4.4 below).

In the scaling limit, we set

$$\begin{align} \displaystyle n = \frac{L y}{\ell}, \nonumber \end{align} \tag{ 61 }$$

under which $P_{\sigma_1\sigma_2}(y) = \frac{L}{\ell} P_{\sigma_1\sigma_2}(n)$. The discrete distribution (44) contains a set of terms of the form

$$\begin{align} \displaystyle \label{dterm} a_{\sigma_{1}\sigma_{2}}(z_{\rho}) z_{\rho}^{-n+1} \quad {\rm~where~}\quad a_{\sigma_{1}\sigma_{2}}(z_{\rho}) = \frac{-J_{\sigma_{1}\sigma_{2}}(z_{\rho})}{[q(x)/(x-z_{\rho})]_{|x=z_{j=\rho}}} \nonumber \end{align} \tag{ 62 }$$

and $z_{\rho}$ is one of the five roots, $z_{\rho} \in \{1, z_+, 1/z_+, z_-, 1/z_-\}$. We now establish their behaviour in the scaling limit.

The easiest case to deal with $z_{\rho}=1$, where the amplitude $a_{\sigma_1\sigma_2}$ that appears in the result for the scaling limit, (5)–(8), is equal to $\lim\nolimits_{L\to\infty}\frac{L}{\ell} a_{\sigma_1\sigma_2}(1)$.

At $z_{\rho}=z_+$ we have the combination

$$\begin{align} \displaystyle \lim_{L \to \infty} \frac{L}{\ell} a_{\sigma_{1}\sigma_{2}}(z_+) \left(1 + \frac{\sqrt{2\phi}}{\sqrt{L}}\right)^{-n+1}, \nonumber \end{align} \tag{ 63 }$$

which, in terms of the continuous coordinate y, becomes

$$\begin{align} \displaystyle \lim_{L \to \infty} \frac{L}{\ell} a_{\sigma_{1}\sigma_{2}}(z_+) \exp \left(- \frac{\sqrt{2\phi L}}{\ell} y \right) = \left[\lim_{L \to \infty} a_{\sigma_{1}\sigma_{2}}(z_+) \sqrt{\frac{L}{2\phi}} \right] \delta(y). \nonumber \end{align} \tag{ 64 }$$

The scaling of $a_{\sigma_{1}\sigma_{2}}(z_+)$ in the large L limit determines whether the delta function actually appears in the $\sigma_1\sigma_2$ sector. In particular, if $a_{\sigma_{1}\sigma_{2}}(z_+)$ decays faster than $1/\sqrt{L}$, we will not get a delta function contribution. The quantity in the square bracket can be identified as $b_{\sigma_1\sigma_2}^{(0)}$ that appears in the probability distribution in the scaling limit. Note that in equations (5)–(8) we dropped superscripts on the amplitudes in the scaling limit where this was unambiguous.

At $z_\rho = 1/z_+$, we find

$$\begin{align} \displaystyle \left[\lim_{L \to \infty} a_{\sigma_{1}\sigma_{2}}(1/z_+) \sqrt{\frac{L}{2\phi}} \right] \delta(\ell-y). \nonumber \end{align} \tag{ 65 }$$

Here, the term in square brackets defines the amplitude $b_{\sigma_1\sigma_2}^{(1)}$.

Turning now to the root $z_{\rho}=z_{-}$, we have

$$\begin{align} \displaystyle \lim_{L \to \infty} \frac{L}{\ell} a_{\sigma_{1}\sigma_{2}}(z_-) \left(1 + \frac{\lambda}{L}\right)^{-n+1} = \left[\lim_{L\to\infty} \frac{L}{\ell} a_{\sigma_{1}\sigma_{2}}(z_-) \right] \exp\left(- \frac{y}{\kappa} \right) \nonumber \end{align} \tag{ 66 }$$

in which we have introduced the lengthscale

$$\begin{align} \displaystyle \kappa = \frac{\ell}{\lambda} = \frac{\sqrt{2} \ell}{\sqrt{(\phi+\theta)(2\phi+\theta)}}. \nonumber \end{align} \tag{ 67 }$$

The square-bracketed term defines the amplitude $c_{\sigma_1\sigma_2}^{(0)}$.

Finally, at $z_{\rho}=1/z_{-}$, we find

$$\begin{align} \displaystyle \lim_{L \to \infty} a_{\sigma_{1}\sigma_{2}}(1/z_-) \left(1 - \frac{\lambda}{L}\right)^{-n+1} = \left[\lim_{L\to\infty} a_{\sigma_{1}\sigma_{2}}(1/z_-) {\rm e}^{\lambda} \right] \exp\left(- \frac{\ell-y}{\kappa} \right), \nonumber \end{align} \tag{ 68 }$$

which furnishes an expression for the amplitude $c_{\sigma_1\sigma_2}^{(1)}$.

It now remains to evaluate the amplitudes. Recall that $J_{\sigma_1\sigma_2}(x)$, defined by (41) is by construction a polynomial of degree $\leqslant 4$. Specifically,

$$\begin{align} \newcommand{\adj}{{{\rm adj}}} \displaystyle J_{\sigma_1\sigma_2}(x) = \hat{T}_4 \sum_j x^2 \adj A_{\sigma_1\sigma_2,j}(x) b_j(x) \nonumber \end{align} \tag{ 69 }$$

where the operator $\hat{T}_4$ discards terms of order x⁵ and higher in a power series in x. This we may write as

$$\begin{align} \newcommand{\tum}{{0}} \displaystyle J_{\sigma_{1}\sigma_{2}}(x) = \tilde{A}_{\sigma_{1}\sigma_{2},1}^{\prime} P_{++}(1) + \tilde{A}_{\sigma_{1}\sigma_{2},2}^{\prime}P_{+-}(1) + \tilde{A}_{\sigma_{1}\sigma_{2},5}^{\prime}P_{+\tum}(1), \nonumber \end{align} \tag{ 70 }$$

where $\newcommand{\adj}{{{\rm adj}}} \tilde{A}_{\sigma_{1}\sigma_{2}, j}^{\prime}(x) = \hat{T}_4 x^2 \adj A_{\sigma_1\sigma_2, j}(x)$.

In the $++$ sector, the adjugate elements exhibit the symmetries

$$\begin{align} \displaystyle \tilde{A}_{++,1}(x) = \tilde{A}_{++,1}(1/x) \nonumber \end{align} \tag{ 71 }$$

$$\begin{align} \displaystyle \tilde{A}_{++,2}(x) = \tilde{A}_{++,3}(1/x) \nonumber \end{align} \tag{ 72 }$$

$$\begin{align} \displaystyle \tilde{A}_{++,5}(x) = \tilde{A}_{++,4}(1/x) \nonumber \end{align} \tag{ 73 }$$

as can be verified by inspection of the explicit expressions presented in the supplementary material. Using these symmetries in (49), one can show that

$$\begin{align} \displaystyle J_{++}(1/z_{\rho}) = z_{\rho}^{1-L}J_{++}(z_\rho) \label{++ symm} \nonumber \end{align} \tag{ 74 }$$

at each of the roots $z_\rho$. The same symmetry also applies in the $\newcommand{\tum}{{0}} \tum\tum$ sector, namely $\newcommand{\tum}{{0}} J_{\tum\tum}(1/z_{\rho}) = z_{\rho}^{1-L}J_{\tum\tum}(z_\rho)$.

Meanwhile, the denominator $[q(x)/(x-z_{\rho})]_{|x=z_{\rho}}$ that appears in (62), has limiting expressions that are symmetric in $z_{\rho} \to 1/z_{\rho}$. These expressions are

$$\begin{align} \displaystyle {[q(x)/(x-1)]}|_{x=1} \sim \frac{\phi (\theta +\phi) (\theta +2 \phi)}{L^3} \label{q1} \nonumber \end{align} \tag{ 75 }$$

$$\begin{align} \displaystyle {[q(x)/(x-z_{+})]}|_{x=z_{+}} \sim {[q(x)/(x-1/z_{+})]}|_{x=1/z_{+}} \sim \frac{8 \phi^{2}}{L^{2}} \label{qzplus} \nonumber \end{align} \tag{ 76 }$$

$$\begin{align} \displaystyle {[q(x)/(x-z_{-})]}|_{x=z_{-}} \sim {[q(x)/(x-1/z_{-})]}|_{x=1/z_{-}} \sim \frac{2 \phi (\theta +\phi) (\theta +2 \phi)}{L^3}.\label{qzminus} \nonumber \end{align} \tag{ 77 }$$

The consequence of these symmetries is that the amplitudes $b_{++}^{(0)}=b_{++}^{(1)} \equiv b_{++}$, $c_{++}^{(0)}=c_{++}^{(1)} \equiv c_{++}$, and similarly $\newcommand{\tum}{{0}} b_{\tum\tum}^{(0)}=b_{\tum\tum}^{(1)} \equiv b_{\tum\tum}$, $\newcommand{\tum}{{0}} c_{\tum\tum}^{(0)}=c_{\tum\tum}^{(1)} \equiv c_{\tum\tum}$.

The remaining ingredient in the amplitudes is the leading large-L behaviour of the truncated adjugate elements $\tilde{A}_{\sigma_{1}\sigma_{2}, j}^{\prime}(x)$ in the scaling limit. In the $++$ sector, these coincide with the expressions set out in table 1. The expressions that are required in the $+-$, $\newcommand{\tum}{{0}} +\tum$ and $\newcommand{\tum}{{0}} \tum\tum$ sectors are provided in tables 2–4.

Table 2. Adjugate elements in the scaling limit for $J_{+-}$.

	$x=1$	$x=z_{+}$	$x=1/z_{-}$	$x=z_{-}$
$\tilde{A}_{+-, 1}$	$-\frac{\theta ^2 \phi ^2 \zeta}{4 L^5}$	$-\frac{\theta ^2 \phi ^2}{L^4}$	$-\frac{\theta ^2 \phi \zeta \left(2\lambda+\zeta \right)}{4 L^5}$	$-\frac{\theta^{2}\phi \zeta(\zeta - 2\lambda)}{4L^{5}}$
$\tilde{A}_{+-, 2}^{\prime}$	$-\frac{\theta ^2 \phi ^2 \zeta}{4 L^5}$	$\frac{\theta ^3 \phi ^{3/2}}{L^{9/2}}$	$\frac{\theta \phi ^2 \zeta \left(4\lambda+\zeta + 2\eta \right)}{4 L^5}$	$\frac{\theta \phi ^2 \zeta \left(-4 \lambda+\zeta + 2\eta \right)}{4 L^5}$
$\tilde{A}_{+-, 5}^{\prime}$	$-\frac{\theta ^2 \phi ^2 \zeta}{4 L^5}$	$\frac{-\sqrt{2} \theta ^3 \phi ^{3/2}}{4 L^{9/2}}$	$-\frac{\theta ^2 \phi \zeta \left(6\lambda+2\zeta + 2\eta \right)}{8 L^5}$	$-\frac{\theta ^2 \phi \zeta \left(-6\lambda+2\zeta + 2\eta \right)}{8 L^5}$

Table 3. Adjugate elements in the scaling limit in $\newcommand{\tum}{{0}} J_{+\tum}$.

	$x=1$	$x=1/z_{-}$	$x=z_{-}$
$\newcommand{\tum}{{0}} \tilde{A}_{+\tum, 1}$	$-\frac{\theta \phi ^3\zeta}{2 L^5}$	$\frac{\theta ^2 \phi \zeta \left(2\lambda+2\zeta \right)}{4 L^5}$	$\frac{\theta ^2 \phi \zeta \left(-2\lambda+2 \zeta \right)}{4 L^5}$
$\newcommand{\tum}{{0}} \tilde{A}_{+\tum, 2}^{\prime}$	$-\frac{\theta \phi ^3 \zeta}{2 L^5}$	$-\frac{\theta \phi ^2 \zeta \left(6\lambda+2\zeta +2\eta \right)}{4 L^5}$	$-\frac{\theta \phi ^2 \zeta \left(-6\lambda+2\zeta + 2\eta \right)}{4 L^5}$
$\newcommand{\tum}{{0}} \tilde{A}_{+\tum, 5}^{\prime}$	$-\frac{\theta \phi ^3 \zeta}{2 L^5}$	$\frac{\theta ^2 \phi \zeta \left(4\lambda+2\zeta + \eta \right)}{4 L^5}$	$\frac{\theta ^2 \phi \zeta \left(-4\lambda+2\zeta + \eta \right)}{4 L^5}$

Table 4. Adjugate elements in the scaling limit for $\newcommand{\tum}{{0}} J_{\tum \tum}$.

	$x=1$	$x=1/z_{-}$
$\newcommand{\tum}{{0}} \tilde{A}^{\prime}_{\tum\tum, 1}$	$-\frac{\phi ^4 \zeta}{L^5}$	$\frac{\theta \phi ^2 \eta (\zeta + 2 \phi)}{2 L^5}$
$\newcommand{\tum}{{0}} \tilde{A}^{\prime}_{\tum\tum, 2}$	$-\frac{\phi ^4 \zeta}{L^5}$	$-\frac{\phi ^3 \zeta\left(2\lambda+\zeta \right)}{L^5}$
$\newcommand{\tum}{{0}} \tilde{A}^{\prime}_{\tum\tum, 5}$	$-\frac{\phi ^4 \zeta}{L^5}$	$\frac{\theta \phi ^2 \zeta \left(2\lambda+2\zeta \right)}{2 L^5}$

Putting this all together, we obtain explicit expressions for the amplitudes that appear in the scaling limit of the stationary probability distribution, equations (5)–(8). The amplitudes that remain finite in the $L\to\infty$ limit are

$$\begin{align} \newcommand{\tum}{{0}} \displaystyle a_{++} = \frac{\zeta \theta ^2 P_{+\tum}(1) \left({\rm e}^{\lambda} (\eta +\lambda)-\eta +\lambda \right)}{4 \eta \ell L \left(\zeta \left({\rm e}^{\lambda}-1\right)+2 \left({\rm e}^{\lambda}+1\right) \lambda \right)} \nonumber \end{align} \tag{ 78 }$$

$$\begin{align} \newcommand{\tum}{{0}} \displaystyle b_{++} = \frac{\theta ^2 P_{+\tum}(1)}{4 L \phi} \nonumber \end{align} \tag{ 79 }$$

$$\begin{align} \newcommand{\tum}{{0}} \displaystyle c_{++} = \frac{\zeta \theta ^3 {\rm e}^{\lambda} P_{+\tum}(1)}{8 \ell L \phi \left({\rm e}^{\lambda} (\eta +\lambda)+\eta -\lambda \right)} \nonumber \end{align} \tag{ 80 }$$

$$\begin{align} \newcommand{\tum}{{0}} \displaystyle a_{+-} = \frac{\zeta \theta ^2 P_{+\tum}(1)\left({\rm e}^{\lambda} (\eta +\lambda)-\eta +\lambda \right)}{4 \eta \ell L \left(\zeta \left({\rm e}^{\lambda}-1\right)+2 \left({\rm e}^{\lambda}+1\right) \lambda \right)} \nonumber \end{align} \tag{ 81 }$$

$$\begin{align} \newcommand{\tum}{{0}} \displaystyle c_{+-}^{(0)} = \frac{\zeta \theta ^3 {\rm e}^{\lambda} P_{+\tum}(1)}{8 \ell L \left({\rm e}^{\lambda} \left(2 \theta ^2+3 \theta (\lambda +2 \phi)+4 \phi (\lambda +\phi)\right)+\theta \lambda \right)} \nonumber \end{align} \tag{ 82 }$$

$$\begin{align} \newcommand{\tum}{{0}} \displaystyle c_{+-}^{(1)} = -\frac{\zeta \theta ^2 {\rm e}^{\lambda} P_{+\tum}(1) (\eta +\lambda)}{4 \eta \ell L \left(\zeta \left({\rm e}^{\lambda}-1\right)+2 \left({\rm e}^{\lambda}+1\right) \lambda \right)} \nonumber \end{align} \tag{ 83 }$$

$$\begin{align} \newcommand{\tum}{{0}} \displaystyle a_{+\tum} = \frac{\zeta \theta P_{+\tum}(1) \phi \left({\rm e}^{\lambda} (\eta +\lambda)-\eta +\lambda \right)}{2 \eta \ell L \left(\zeta \left({\rm e}^{\lambda}-1\right)+2 \left({\rm e}^{\lambda}+1\right) \lambda \right)} \nonumber \end{align} \tag{ 84 }$$

$$\begin{align} \newcommand{\tum}{{0}} \displaystyle c_{+\tum}^{(0)} = \frac{\zeta \theta ^2 {\rm e}^{\lambda} P_{+\tum}(1) (2 \lambda -\eta)}{4 \eta \ell L \left(\zeta \left({\rm e}^{\lambda}-1\right)+2 \left({\rm e}^{\lambda}+1\right) \lambda \right)} \nonumber \end{align} \tag{ 85 }$$

$$\begin{align} \newcommand{\tum}{{0}} \displaystyle c_{+\tum}^{(1)} = \frac{\zeta \theta ^2 {\rm e}^{\lambda} P_{+\tum}(1) (\eta +2 \lambda)}{4 \eta l L \left(\zeta \left({\rm e}^{\lambda}-1\right)+2 \left({\rm e}^{\lambda}+1\right) \lambda \right)} \nonumber \end{align} \tag{ 86 }$$

$$\begin{align} \newcommand{\tum}{{0}} \displaystyle a_{\tum \tum} = \frac{\zeta P_{+\tum}(1) \phi ^2 \left({\rm e}^{\lambda} (\eta +\lambda)-\eta +\lambda \right)}{\eta \ell L \left(\zeta \left({\rm e}^{\lambda}-1\right)+2 \left({\rm e}^{\lambda}+1\right) \lambda \right)} \nonumber \end{align} \tag{ 87 }$$

$$\begin{align} \newcommand{\tum}{{0}} \displaystyle c_{\tum \tum} = \frac{\zeta \theta {\rm e}^{\lambda} P_{+\tum}(1) \phi}{2 \ell L \left({\rm e}^{\lambda} (\eta +\lambda)+\eta -\lambda \right)}. \nonumber \end{align} \tag{ 88 }$$

All other amplitudes are zero.

The w coefficients are more straightforward to obtain, since the Kronecker delta symbols $\delta_{n, 1}$ and $\delta_{n, L-1}$ in (32) turn into Dirac delta functions $\delta(y)$ and $\delta(\ell-y)$, respectively, with their amplitudes unchanged. These amplitudes, $w_{\sigma_1\sigma_2}^{(0)}$ and $w_{\sigma_1\sigma_2}^{(1)}$ are found to be

$$\begin{align} \newcommand{\tum}{{0}} \displaystyle w_{+-}^{(0)} = \frac{1}{2}(\alpha+2)\beta P_{+-}(1) + \beta^{2}P_{+-}(1) + \beta^{2}P_{+\tum}(1)/2 \nonumber \end{align} \tag{ 89 }$$

$$\begin{align} \newcommand{\tum}{{0}} \displaystyle \sim \frac{P_{+-} (1)\theta}{L} = \frac{\theta ^2 P_{+\tum}(1) \left({\rm e}^{\lambda} (\zeta +\lambda)-\zeta +\lambda \right)}{L \phi \left(\zeta \left({\rm e}^{\lambda}-1\right)+2 \left({\rm e}^{\lambda}+1\right) \lambda \right)} \nonumber \end{align} \tag{ 90 }$$

$$\begin{align} \newcommand{\tum}{{0}} \displaystyle w_{+\tum}^{(0)} = \alpha\beta P_{+-}(1) + \alpha \beta P_{+\tum}(1) + \beta P_{+\tum}(1) \sim \frac{P_{+\tum}(0) \theta}{L} \nonumber \end{align} \tag{ 91 }$$

$$\begin{align} \newcommand{\tum}{{0}} \displaystyle w_{\tum \tum}^{(0)} = \alpha^{2} P_{+-}(1) +(\alpha + \alpha^{2}) P_{+\tum}(1) \sim \frac{P_{+\tum}(1) \phi}{L} \nonumber \end{align} \tag{ 92 }$$

$$\begin{align} \newcommand{\tum}{{0}} \displaystyle w_{\tum \tum}^{(1)} = w_{\tum \tum}^{(0)}. \nonumber \end{align} \tag{ 93 }$$

Again the other w amplitudes are all zero.

Note that although all the amplitudes have the superficial appearance of a $1/L$ decay, this is in fact cancelled by the remaining constant, $\newcommand{\tum}{{0}} P_{+\tum}(1)$, which scales as L (as will be determined below by normalising the distribution). There is now one final remaining constant, $\newcommand{\tum}{{0}} P_{+\tum}(1)$, which is fixed by normalisation.

Rather than impose normalisation on the whole probability distribution, it is sufficient (and more straightforward) to impose it on a single velocity sector in order to determine $\newcommand{\tum}{{0}} P_{+\tum}(1)$. The relative weight of each sector can be calculated straightforwardly because transitions between sectors occur at rates that are decoupled from the hopping dynamics i.e. the transitions between sectors are independent of the particle separation n. Moreover, each particle enters a velocity state independently of the other. Consequently, if we define the marginal probability distribution

$$\begin{align} \displaystyle P_{\sigma_1} = \sum_{\sigma_2} \int_0^{\ell} {\rm d} y\, P_{\sigma_1\sigma_2}(y) \nonumber \end{align} \tag{ 94 }$$

then we have for the probability of being in the velocity sector $\sigma_1\sigma_2$ that

$$\begin{align} \displaystyle P_{\sigma_1\sigma_2} = \int_0^{\ell} {\rm d} y\, P_{\sigma_1\sigma_2}(y) = P_{\sigma_1} P_{\sigma_2}. \nonumber \end{align} \tag{ 95 }$$

The master equation for the single particle velocity distribution reads

$$\begin{align} \newcommand{\tum}{{0}} \displaystyle \frac{\partial P_{+}}{\partial t} = -\alpha P_{+} + \frac{\beta}{2} P_{\tum} \nonumber \end{align} \tag{ 96 }$$

$$\begin{align} \newcommand{\tum}{{0}} \displaystyle \frac{\partial P_{-}}{\partial t} = -\alpha P_{-} + \frac{\beta}{2} P_{\tum} \nonumber \end{align} \tag{ 97 }$$

$$\begin{align} \newcommand{\tum}{{0}} \displaystyle \frac{\partial P_{\tum}}{\partial t} = \alpha [P_{+} + P_{-}] - \beta P_{\tum}. \nonumber \end{align} \tag{ 98 }$$

In the steady state, we have $P_{+} = P_{-}$ by symmetry and consequently

$$\begin{align} \newcommand{\tum}{{0}} \displaystyle P_{\tum} = \frac{2\alpha}{\beta} P_{+}. \nonumber \end{align} \tag{ 99 }$$

Using this result and the fact that $\newcommand{\tum}{{0}} P_++P_\tum+P_-=1$, we find

$$\begin{align} \displaystyle P_{+} = P_{-} = \frac{1}{2(1+\alpha/\beta)}. \nonumber \end{align} \tag{ 100 }$$

Insisting now that $\int_0^{\ell} {\rm d}\ell P_{++}(y) = P_{+}^2$, we find that

$$\begin{align} \newcommand{\tum}{{0}} \displaystyle P_{+\tum}(1) =L\left[(\theta +\phi)^2 \left(\zeta \left(\frac{\sqrt{2} \theta \left({\rm e}^{\lambda}-1\right)}{\phi \sqrt{\zeta \eta} \left(\eta +{\rm e}^{\lambda} (\theta +\lambda +\phi)-\lambda \right)}+\frac{{\rm e}^{\lambda} (\eta +\lambda)-\theta +\lambda -\phi}{\eta \left(\zeta \left({\rm e}^{\lambda}-1\right)+2 \left({\rm e}^{\lambda}+1\right) \lambda \right)}\right)+\frac{2}{\phi}\right)\right]^{-1}, \nonumber \end{align} \tag{ 101 }$$

which completes our derivation of equations (5)–(8).

We can directly simulate the scaling limit by having particles move ballistically at speed v and undergoing tumbling and untumbling events at times drawn from an exponential distribution with means $1/\alpha$ and $1/\beta$ respectively. In figures 3 and 4 we compare the distributions (in the form of effective potentials) obtained from this simulation with our analytical calculation. Once again, we find complete agreement. As discussed in section 1.2, and as seen explicitly above, $\xi_{+}$ collapses to a delta function in this limit. Nevertheless, the second lengthscale, $\xi_{-}$, which is induced by the finite tumbling time, remains physically relevant in the scaling limit.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As a further check, we may consider the limit where the exit rate from tumbling $\tilde{\beta} \rightarrow \infty$. In this limit tumbling is instantaneous, and we recover the probability distribution in the scaling limit of the model studied in [11]

$$\begin{align} \displaystyle P_{++}(y) = \frac{\phi}{4\ell(4+\phi)} + \frac{\delta (y) + \delta (\ell - y)}{2(4+ \phi)} = \frac{\tilde\alpha + 2v [\delta(y) + \delta(\ell-y)]}{4(4v+\tilde\alpha\ell)} \label{++ beta inf} \nonumber \end{align} \tag{ 102 }$$

$$\begin{align} \displaystyle P_{+-}(y) = \frac{\phi}{4\ell (4+\phi)} + \frac{\delta (y)}{(4+\phi)} = \frac{\tilde\alpha + 4v \delta(y)}{4(4v+\tilde\alpha\ell)} \nonumber \end{align} \tag{ 103 }$$

$$\begin{align} \newcommand{\tum}{{0}} \displaystyle P_{+\tum}(y) = P_{\tum \tum}(y)= 0. \nonumber \end{align} \tag{ 104 }$$

Moreover, it is instructive to note exactly how this limit is recovered.

As expected, all contributions from states with a tumbling particle vanish in this limit. There are no contributions from $c_{+-}^{(0)}{\rm e}^{-y/\kappa}$ and $c_{+-}^{(1)}{\rm e}^{-(\ell-y)/\kappa}$ as the exponentials vanish. Therefore the only terms that contribute from $(+-)$ are

$$\begin{align} \displaystyle a_{+-} \sim \frac{\phi}{4\ell (4+\phi)} \;\; {\rm~and,~} \nonumber \end{align} \tag{ 105 }$$

$$\begin{align} \displaystyle w_{+-}^{(0)} \sim \frac{1}{4 (4+\phi)} \;\cdot \nonumber \end{align} \tag{ 106 }$$

However, all of the terms in $(++)$ (and, equivalently, its symmetric counterpart $(--)$) do contribute to the probability in this limit. Specifically, the constant

$$\begin{align} \displaystyle a_{++} \sim \frac{\phi}{4\ell (4+\phi)}, \nonumber \end{align} \tag{ 107 }$$

and

$$\begin{align} \displaystyle b_{++}[\delta(y) + \delta(\ell -y)] \sim \frac{\sqrt{2}}{4 (4+\phi)} [\delta(y) + \delta(\ell -y)] \;\;\; {\rm~and~} \nonumber \end{align} \tag{ 108 }$$

$$\begin{align} \displaystyle c_{++}[{\rm e}^{- y/\kappa} + {\rm e}^{-(\ell-y)/\kappa}] \sim \left(\frac{1}{2(4+\phi)} - \frac{\sqrt{2}}{4(4+\phi)} \right) [\delta(y) + \delta(\ell -y)]. \nonumber \end{align} \tag{ 109 }$$

Thus we see that not all of the probability in the delta functions in (102) comes from the delta-function term $b_{++}$, but that there is also a contribution from the originally finite exponential piece multiplying $c_{++}$. In other words, when the tumbling time is short (but not zero), very small inter-particle separations are generated with a high probability as a consequence of the short distance moved by a particle following a collision while the other one is tumbling. The effect of this is seen in simulations: when β is set very large but not strictly infinite there is a significant fraction of the probability for configurations at very marginal but non-zero separations. Only when β is set strictly infinite does this probability moves into the delta-function terms.