### Appendix A: variation of distribution with model parameters

To better understand the range of behavior possible for the model, we relax the requirement that parameters must be fixed to measured values. As described above, the density in the interior is constant by symmetry. It formally diverges at the boundary since the current of ants incident on the \(x=L/2\) boundary during a given time step will all be placed at \(x=L/2\) for the next time step: a finite number of members occupying a single point. The value of the basic distance scale \(\ell\), over which *n*(*x*) varies, changes as we modify parameters: larger mean-square impulses \(\sigma _L^2\) and a smaller drag force coefficient \(1/\tau _D\) increase the region that a given ant explores, extending the range over which boundaries influence the distribution. But the qualitative behavior of the *n*(*x*) does not change dramatically.

The distinctive shape of \(P(v_x)\) with round shoulder peaks coming from the interior and a sharp central peak from the boundaries, is sensitive to parameter values. In Fig. 5 we show the distribution for various parameter choices. In the upper panels we vary \(\tau _D\) and \(v_0\), which control the deterministic force, and in the lower panels vary \(\sigma _L\) and \(\sigma _T\), the size of the random impulses.

The shoulder peaks are centered on \(\pm v_0\) and have a width \(v_\infty\). Increasing the driving force \(v_0\) pushes them outward without changing their height. Upon decreasing \(v_0\), the shoulders can merge with the central peak, so the width of the latter increases significantly. If the drag coefficient \(1/\tau _D\) increases then accumulated random impulses cannot achieve large speeds so \(v_\infty\) is decreased. The result is the shoulder peaks become narrow, tall and more pronounced. On the other hand, members become more narrowly localized in space so \(\ell\) decreases and the boundaries affect a smaller region. The height of the central peak slightly decreases as a result. If \(1/\tau _D\) decreases and \(v_\infty\) gets large then the shoulders become wider and shorter, merging into a smooth distribution.

With the deterministic force fixed we can vary the size of the root-mean-square random impulses \(\sigma _L\) and \(\sigma _T\). Variation of \(\sigma _L\) has a similar effect to changing \(\tau _D\), as is clear from the two plots on the left of Fig. 5. This is because the distribution results from a competition between the deterministic force, which scales as \(1/\tau _D\), and the random impulses that scale as \(\sigma _L^2\). So, increasing one of these is similar to decreasing the other. Recall, in Eq. 19, the anisotropy of the impulses produced only logarithmic corrections to the distribution in the interior. As a result, when we change \(\sigma _T\) the shoulder features look the same, as seen in the lower right of Fig. 5. But the central peak is more sensitive to \(\sigma _T\). After a model ant hits the boundary its body is aligned parallel with it (because the normal component of its velocity is reset to zero). If \(\sigma _T\) is small, it cannot change direction effectively so it moves along the edge and remains in contact with the boundary. The height of the central peak, which indicates the weight of the boundary population, increases as a result.

### Appendix B: Equilibrium distribution for a 1D arena

Since the rate equation for the bounded 2D model is more difficult to solve, we consider a simpler 1D model in this section. The solution for the 1D equilibrium distribution reproduces the key properties of the density *n*(*x*) and velocity distribution \(P(v_x)\) found above in the 2D simulations and experiment. The 1D model can be viewed as the limiting case of extremely anisotropic impulses for the 2D model. We set \(\sigma _T=0\) so an ant initially oriented along the *x* axis remains so.

The arena is one-dimensional with length \(R=L/\ell\) centered on the origin \(x=0\), with parameters for speed \(v_\infty\) and distance \(\ell\) defined in Eq. 16. The position and velocity are *x* and *v* have dimensionless counterparts \(r=x/\ell\) and \(u=v/v_\infty\). The force *F* and its dimensionless version *f* are

$$\begin{aligned} F=-\frac{v}{\tau _D}+\frac{v_0\textrm{sign}(v)}{\tau _D},\;\;\;\;f=-u+u_0\textrm{sign}(u) \end{aligned}$$

(23)

and the rate equation that determines the equilibrium distribution \(\Pi (r,u)\) is

$$\begin{aligned} 0=\frac{\partial }{\partial u}\bigg (-f\Pi +\frac{\partial }{\partial u}\Pi \bigg )-u\frac{\partial }{\partial r}\Pi . \end{aligned}$$

(24)

Following the protocol for the simulations described above, we know that the boundary distribution \(\Pi (R/2,u)=\Pi (-R/2,-u)\) should be similar to \(\delta (u)\) because all members that reach a boundary are restarted with zero velocity. However, at a boundary, the distribution is not an even function of *u*.

To solve Eq. 24 it is helpful first to find eigenfunctions \(g_\lambda (u)\) of the velocity operator appearing in it, in the sense that

$$\begin{aligned} \bigg (-f g_\lambda (u)+g’_\lambda (u)\bigg )’=\lambda u g_\lambda (u),\;\;\textrm{with}\;g_\lambda (0)=1 \end{aligned}$$

(25)

for a given eigenvalue \(\lambda\). (The primes denote derivatives with respect to *u*.) Having found such functions, we can write the general solution for the equilibrium distribution as a sum over \(g_\lambda (u)\) functions weighted by coefficients \(c_\lambda (r)\) that vary exponentially with position,

$$\begin{aligned} \Pi (r,u)=\sum _\lambda c_\lambda (r)g_\lambda (u)\;\;\textrm{where}\; c_\lambda (r)=c_\lambda (R/2)\textrm{e}^{\lambda [r-R/2]}. \end{aligned}$$

(26)

Also, Eq. 25 implies \(g_{-\lambda }(u)=g_\lambda (-u)\) and the equilibrium distribution must satisfy reflection symmetry so \(\Pi (r,u)=\Pi (-r,-u)\). This means we can restrict the sum to non-negative eigenvalues and write

$$\begin{aligned} \Pi (r,u)=\sum _{\lambda \ge 0} c_\lambda (R/2)\bigg (g_\lambda (u)\textrm{e}^{\lambda [r-R/2] }+g_\lambda (-u)\textrm{e}^{-\lambda [r+R/2]}\bigg ). \end{aligned}$$

(27)

The \(\lambda =0\) term is position-independent: it gives the equilibrium distribution for an unbounded arena. The \(\lambda >0\) terms vary with position: the first term on the right of Eq. 27 dominates near the right boundary but is negligible at the left boundary. The second term has the opposite property.

The eigenvalue equation, Eq. 25, for \(\lambda =0\) requires that

$$\begin{aligned} -f(u)g_0(u)+g_0′(u)=-f(\infty )g_0(\infty )+g’_0(\infty )=0, \end{aligned}$$

(28)

where the last equality is required for a (normalizable) physical distribution. It is solved by

$$\begin{aligned} g_0(u)=\exp (\chi _0[u])\;\;\textrm{with}\;\;\chi _0(u)=\int _0^u du’ f(u’)=-\frac{u^2}{2}+|u|u_0. \end{aligned}$$

(29)

Clearly, \(g_0(u)\) is the 1D analogue of the equilibrium distribution found above, in Eq. 19, for the 2D unbounded arena. According to Eq. 27, the distribution approaches \(g_0(u)\) as one moves further from a boundary.

The eigenfunction \(g_\lambda (u)\) for \(\lambda \ne 0\) can be expressed exactly as a Hermite function, with an order that increases with \(\lambda\), multiplied by a Gaussian factor. It oscillates increasingly rapidly as \(\lambda\) increases, so a \(g_\lambda (u)\) function is analogous to a Fourier component. Near the boundary, many \(g_\lambda (u)\) contribute and their sum gives a velocity distribution that is sharply peaked near \(u=0\).

To make it easier to sum over \(g_\lambda (u)\) with \(\lambda >0\) we approximate these functions. Consider Eq. 25 when \(\lambda>>1\) and the terms expected to dominate it satisfy \(g”_\lambda (u)-\lambda u g_\lambda (u)=0\). This is solved by an Airy function of the first kind: \(g_\lambda (u)\propto \textrm{Ai}(\lambda ^{1/3}u)\). (An Airy function of the second kind also solves this equation but would give a physically incorrect distribution.) We substitute a trial form \(g_\lambda (u)=\textrm{Ai}(\lambda ^{1/3}u)\exp (\chi _1[u])\), with \(\chi _1(u)\) an unknown function, into Eq. 25 and obtain

$$\begin{aligned} 0=\lambda ^{1/3}\textrm{Ai}'(\lambda ^{1/3}u)\bigg (-f(u)+2\chi _1′(u)\bigg )+\textrm{Ai}(\lambda ^{1/3}u)\bigg (-f'(u)-f(u)\chi _1′(u)+\chi _1”(u)+[\chi _1′(u)]^2\bigg ). \end{aligned}$$

(30)

For sufficiently large \(\lambda\) the first expression in parentheses, which is multiplied by \(\lambda ^{1/3}\), must vanish. It does if \(\chi _1(u)=\chi _0(u)/2\). The result

$$\begin{aligned} g_\lambda (u)=\exp (\chi _0[u]/2)\frac{\textrm{Ai}(\lambda ^{1/3}u)}{\textrm{Ai}(0)} \end{aligned}$$

(31)

agrees reasonably well with the true eigenfunction for \(\lambda \approx 1\) and becomes exact in the limit of large \(\lambda\). It is convenient because it is well behaved for any eigenvalue so we can treat \(\lambda >0\) as a continuous parameter. Note that \(\textrm{Ai}(0)=(3^{2/3}\Gamma [2/3])^{-1}\approx 0.36\) and the Airy function is positive and exponentially decaying when its argument is positive but oscillates in sign for negative argument.

We replace the sum over non-zero \(\lambda\) with an integral using an interpolating function \(c(\lambda )\equiv c_\lambda (R/2)\). A minimum value of \(\lambda =1\) will be imposed on the continuous spectrum–the precise minimum value will not affect important results below. The distribution becomes:

$$\begin{aligned} \Pi (r,u)=c_0\textrm{e}^{\chi _0(u)}+ \frac{\textrm{e}^{\chi _0[u]/2}}{\textrm{Ai}(0)}\int _{1}^{\Lambda } d\lambda c(\lambda )\bigg ( \textrm{Ai}(u\lambda ^{1/3})\textrm{e}^{\lambda [r-R/2]}+\textrm{Ai}(-u\lambda ^{1/3})\textrm{e}^{-\lambda [r+R/2]}\bigg ). \end{aligned}$$

(32)

The cutoff of the \(\lambda\) integral can be chosen according to the desired position resolution \(\Delta r\) and velocity resolution \(\Delta u\) of the distribution. The terms with \(\lambda \gg 1/\Delta r\) are zero at any measurable distance from the boundary, so they can be ignored. Also, since \(\textrm{Ai}(\pm \infty )=0\), terms with \(\lambda>>1/\Delta u\) contribute nothing to the distribution. For our, qualitative, purpose we set \(\Lambda =\infty\) whenever the integral converges. Also, if the \(\lambda =0\) term does not dominate the sum then we can extend the lower integration limit to \(\lambda =0\) (it does not matter if we use Eq. 31 to poorly approximate \(g_0(u)\) in this case).

The \(c(\lambda )\) coefficients should be determined by matching the \(r=R/2\) distribution

$$\begin{aligned} \Pi (R/2,u)=\frac{\textrm{e}^{\chi _0(u)/2}}{\textrm{Ai}(0)} \int _{0}^{\Lambda } d\lambda c(\lambda )\textrm{Ai}(u\lambda ^{1/3}),\;\;\;\Pi (R/2,0)=\int _0^\Lambda d\lambda c(\lambda ) \end{aligned}$$

(33)

to a known boundary value. Current conservation at the boundary determines the rate that members are replaced at \(r=R/2\) with zero velocity. But as a rough approximation to this procedure, and in analogy with the corresponding situation for a Fourier transform, we assume that \(u=0\) members at \(r=R/2\) are uniformly distributed among all \(\lambda\) components. This means \(c(\lambda )=c_1\approx c_0\) with \(c_1\) a constant.

The border distribution is now

$$\begin{aligned} \Pi (R/2,u)\approx \frac{3c_1}{\textrm{Ai}(0)u^3}\bigg (\textrm{Ai}(0)-\textrm{Ai}(\Lambda ^{1/3}u)+u\Lambda ^{1/3}\textrm{Ai}'(u\Lambda ^{1/3})\bigg ) \end{aligned}$$

(34)

where we used the fact that the integral is strongly peaked about \(u=0\) to set \(\chi _0(u)\approx \chi _0(0)=0\) and then used the defining differential equation of the Airy function to do the \(\lambda\) integral. At \(u=0\) the function has a value \(\Pi (R/2,0)\approx \Lambda\) and is peaked at the slightly negative velocity of \(u=u_{\max }\approx -1.47\Lambda ^{-1/3}\) with \(\Pi (R/2,u_{\max })\approx 1.46\Lambda\). For positive velocities it drops to zero on a scale of \(u_{\max }\) while for negative velocities it oscillates on this same scale, rapidly averaging to zero. In Fig. 6 we plot it using a modest value \(\Lambda ^{1/3}=20\) so the qualitative behavior can be clearly seen.

The overall velocity distribution is similarly obtained

$$\begin{aligned} P(u)=\int _{-R/2}^{R/2} dr \Pi (r,u)=c_0R\textrm{e}^{\chi _0(u)}+c_1\textrm{e}^{\chi _0(u)/2}B(u) \end{aligned}$$

(35)

where

$$\begin{aligned} B(u)=\frac{3}{\textrm{Ai}(0)}\int _{1}^{\Lambda ^{1/3}} \frac{dz}{z} \bigg ( \textrm{Ai}(zu)+\textrm{Ai}(-zu)\bigg ). \end{aligned}$$

(36)

The function *B*(*u*) comes from the population near the boundaries and determines the velocity distribution at small *u*. It has a \(u=0\) peak of height \(B(0)=2\ln \Lambda\) and width (FWHM) of approximately \(4/\ln \Lambda\). While *B*(*u*) depends on \(\Lambda\), this dependence is weak. In Fig. 6 we plot *B*(*u*) and *P*(*u*) with \(\Lambda =10^6\) and \(u_0=1.5\), \(R=4\) and \(c_0=c_1\). (The vertical axis for plots of *P*(*u*) and *n*(*r*) is arbitrary in this figure: the amplitude of these distributions should be fixed by normalization.) The cusp in *P*(*u*) at \(u=0\) comes from that of the *B*(*u*) function while the shoulders, coming from \(g_0(u)\), correspond to \(u=\pm u_0\).

To obtain the density

$$\begin{aligned} n(r)=\int _{-\infty }^{\infty } du \Pi (r,u)\;\;\textrm{with}\;n(0)=c_0\int _{-\infty }^{\infty } du \textrm{e}^{\chi _o(u)} \end{aligned}$$

(37)

we integrate Eq. 32 numerically. The \(\lambda\) integral converges so we can take \(\Lambda \rightarrow \infty\). The result for \(n(r)-n(0)\) is plotted in Fig. 6 versus the distance from either boundary.

Insight can be gained by considering the case of no driving force, i.e. setting \(u_0=0\), so that analytic results for *n*(*r*) are possible. Integrating Eq. 32 over *u*, using the large \(\lambda\) approximation of \(2\sqrt{\pi }\textrm{e}^{2\lambda ^2/3}\textrm{Ai}(\lambda ^{4/3})\approx \lambda ^{-1/3}\), and taking \(\Lambda \rightarrow \infty\) we find

$$\begin{aligned} n(r)-n(0)\approx c_1\int _1^\infty d\lambda \lambda ^{-1/3} \bigg (\textrm{e}^{\lambda (r-R/2)}+\textrm{e}^{-\lambda (r+R/2)}\bigg )\;\;\;\textrm{for}\;u_0=0. \end{aligned}$$

(38)

Writing the positive distance from the left boundary as \(d=R/2+r\) or from the right boundary as \(d=R/2-r\), we have

$$\begin{aligned} n(r)-n(0)\approx \frac{c_1\Gamma (2/3)}{d^{2/3}}\;\; (\textrm{for}\;d<<1)\;\;,\; n(r)-n(0)\approx \frac{c_1\textrm{e}^{-d}}{d}\;\;(\textrm{for}\;d>>1). \end{aligned}$$

(39)

If we had extended the lower limit of the integral to \(\lambda =0\) then the density, for all distances *d*, would be given by the \(d<<1\) expression above (so the lower cutoff in the \(\lambda\) integral affects the long-distance decay of the boundary population). In Fig. 6 we compare the \(d^{-2/3}\) divergence with the numerical calculation of *n*(*r*). To within a constant multiplicative factor, the \(u_0=0\) approximation of the density is similar as the result for finite \(u_0\).

The results above are consistent with what we saw in the experiment and the model simulations. The density *n*(*x*) is large at the boundary because so many slow-moving ants remain near it. The density decreases with distance *d* and ants further than \(\ell\) from the boundary act as though they are in an infinite arena. The sharp \(v=0\) peak in the velocity distribution *P*(*v*) is due to ants at the boundary that are often abruptly stopped. The shoulder features in *P*(*v*) are properties of the distribution in the arena interior: they are local maxima in the distribution centered on \(v=v_0\) with a width of \(v_\infty\).

## Add Comment