{'fraction_non_alphanumeric': 0.029620162760539026, 'fraction_numerical': 0.008814629653992895, 'mean_word_length': 4.535580524344569, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 9, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Lattice Models of Nonequilibrium Bacterial Dynamics

We study a model of self propelled particles exhibiting run and tumble dynamics on lattice. This non-Brownian diffusion is characterised by a random walk with a finite persistence length between changes of direction, and is inspired by the motion of bacteria such as E. coli. By defining a class of models with multiple species of particle and transmutation between species we can recreate such dynamics. These models admit exact analytical results whilst also forming a counterpart to previous continuum models of run and tumble dynamics. We solve the externally driven non-interacting and zero-range versions of the model exactly and utilise a field theoretic approach to derive the continuum fluctuating hydrodynamics for more general interactions. We make contact with prior approaches to run and tumble dynamics off lattice and determine the steady state and linear stability for a class of crowding interactions, where the jump rate decreases as density increases. In addition to its interest from the perspective of nonequilibrium statistical mechanics, this lattice model constitutes and efficient tool to simulate a class of interacting run and tumble models relevant to bacterial motion, so long as certain conditions (that we derive) are met.


Introduction
The last decade has seen a growing number of studies of biological systems conducted by physicists. Much of this work relies on tools not traditionally found in biology. For instance, recent methods of nonlinear optics have made possible observation of biological systems on much smaller scales than was previously possible [1]. In other cases biological systems have helped shed light on questions of fundamental importance in theoretical physics. Studies of bird flocks showed that it may be possible to observe long range order in systems with a continuous symmetry if detailed balance is broken [2] (in equilibrium this is forbidden by the Mermin-Wagner theorem [3]).
Bacteria provide one example of a biological system which is of great interest to the study of nonequilibrium physics. Some species of bacteria, such as Escherichia coli, swim by means of a series of relatively straight runs interspersed with short periods of chaotic motion, known as tumbles, during which their orientation changes at random and they experience little net displacement [4]. This can be seen as a type of non-Brownian diffusion where the steady state probability distribution will not be Boltzmann.
Following experiments to determine their behaviour by Berg and others [5][6][7][8][9][10][11] since the 1970s, much is known about the dynamics and behaviour of individual bacteria. Less, however, is known about their collective behaviour and it is here that statistical mechanics can play a useful role. Schnitzer et al. analysed various strategies bacteria may employ, e.g. changing their speed or tumble rate, and studied the differences these make to the steady state distribution [12]. These analyses, which were based on off-lattice descriptions of the particle dynamics, were more recently extended to allow for various types of interaction between bacteria [13], external force fields [14], hydrodynamic interactions [15], and coupling to logistic population growth [16]. Intriguingly it was shown that a dependence of motility parameters such as swim-speed on density was alone enough to cause phase separation [13], in contrast to detailed-balance systems where many-body effects on particle mobility cancel exactly in Boltzmann steady states.
In this work we construct a lattice model whose non-interacting continuum limit recovers the off lattice equations previously derived for noninteracting run-and-tumble particles [12,14]. We will also address interactions in the form of density dependent motility parameters, and thereby connect also with the off lattice approach of [13]. One motivation for this approach is as follows. Although real bacterial systems are off lattice, both the microscopic run-and-tumble equations, and the diffusive continuum equations found by coarse graining these, are difficult to simulate efficiently, particularly once interactions are included. On lattice the local density of particles is easy to determine and so density-dependent interactions are easy to include at relatively low cost computationally. Lattice simulations are also easier to extend to higher dimensions as we will consider towards the end of this paper. For all these reasons, creation of a robust and accurate lattice representation of bacterial motility is a worthwhile goal, even from a purely phenomenological standpoint, whereby the purpose of a model is to provide a fairly direct explanation for results seen in experiment.
A second motivation is more fundamental. Models of non-equilibrium statistical mechanics can be broadly split into two categories: one phenomenological as just described, the other comprising simple models which allow the study of basic concepts and facilitate a more detailed understanding of the nature of non-equilibrium physics. In the latter category we can think of lattice transport models such as the exclusion [17] or zero-range [18] processes, for which some exact analytical results can be found, as well as methods to characterise fluctuations and large deviations in non-equilibrium states [19]. Our model falls squarely into this category of simple theoretical models, and indeed extends some of these examples; we investigate both a zero-range process and a partially excluding version of our run-and-tumble system. As a microscopic lattice model we can, under certain conditions, calculate exact steady states and understand precisely how changes in the underlying dynamics affect the probability distributions. More generally we can always write an exact master equation for probability of a configuration and utilise a variety of field theoretic representations to derive the large scale fluctuating hydrodynamics, to attempt to map the system to a free energy and to determine the steady state behaviour and dynamic stability.
As already explained, in addition to being an interesting toy model that we can use to understand driven non-equilibrium physics, our model also forms a direct lattice counterpart to some well established phenomenological continuum models of bacterial dynamics. It is unusual for a model to allow exact computations while credibly describing the real behaviour of a physical system, and this is one of the reasons why we consider this model to be of interest.
In section 2 we introduce our model of run-and-tumble bacteria on a lattice. In section 3 we consider those cases in which we can find exact solutions to the steady state behaviour. After first examining in section 3.1 the non-interacting case, which for inhomogeneous and anisotropic jump and tumble rates can still admit nontrivial solutions, in section 3.2 we include a zero-range interaction and derive sufficient conditions on the rates for a factorised steady state to exist.
In section 4 we consider more general interactions, where the jump rates can depend on the full configuration of the system, and derive a continuum fluctuating hydrodynamics using a field theoretic approach. We first set up the theoretical apparatus and re-derive the results for non-interacting systems within this framework; we then calculate the field equations for the more generic interacting system. In section 5 we consider one specific type of interaction and investigate the effect of crowding. We derive a free-energy-like functional from which we can determine the steady state probability distribution and compare the effects of different methods to coarse grain the local density in the interaction terms. Finally we extend our results to higher dimensions more applicable to real systems in section 5.3.  Figure 1. Presentation of the model. Filled circles represent right moving particles while unfilled circles denote left moving particles. Some of the possible transitions are illustrated on the figure.

The Model
To model the finite persistence length in run-and-tumble dynamics on lattice particles jump repeatedly in the same direction, u, with rate d(u) and change direction -tumblewith rate α(u). In 1d, this means particles can be either right-going or left-going. Rightgoing particles jump with rate d + i from site i to site i+1 and tumble with rate α + i . After a tumble, they become left-going particles with probability 1/2. The corresponding rates for left-going particles are called d − i and α − i , see figure 1. Throughout this paper we assume tumbles to be instantaneous. In reality the duration of tumbles is typically of the order of one tenth of the duration of runs [4]. There may be situations, however, where time spent tumbling may increase, and where the finite duration of a tumble may have an effect on the dynamics and steady state of the system. We leave a consideration of this case to future work.

Exactly Solvable Cases
The jump and tumble rates defined in section 2 may, in general, depend on the occupation numbers of any lattices sites. Through this we may introduce interactions to our system and consider the general case of N interacting bacteria. Before we deal with such complex cases in section 4, however, let us look at the interesting limiting cases that allow for exact computation of the steady state. In particular we start by considering first the non-interacting case and investigating the effect of position and direction dependence on the rates in section 3.1, before then turning to zero-range processes in section 3.2.

Non-Interacting Particles
For non-interacting particles, one can always handle the single particle case first and compute the average occupancy ρ ± i of left (−) or right (+) moving particles on site i. As shown below, the steady-state distribution of n particles is then given by a product on each site of a multinomial distribution of parameters ρ ± i . Let us start with the single particle process and call P (k, +) and P (k, −) the probability to find a bacterium at site k going to the right and to the left, respectively.

The master equation reads
3.1.1. Exact Solution of the Steady State For the case of non-interacting particles, where all interesting and non-trivial effects come from heterogeneities or anisotropies in the jump and tumble rates, we can calculate the steady state probability distributions exactly. We start from the continuity equation where J i,i+1 is the net probability flux between sites i and i + 1, that is As we are only after the steady state distribution we take the time derivative on the left hand side of equation (3) to be zero, so all currents J i,i+1 are constant and equal to some J fixed by the boundary conditions. The master equation for the density of left-moving particles at site i reads Using (4) to eliminate the first term on the right hand side of equation (5), which is the only term to depend on site i + 1, we can establish a relationship between the number of left and right moving particles on site i at steady state: Equation (6), along with equation (4), leads to the recursion relation for right moving particles which can be solved to yield The probability to find a particle at any position and in either state can then be calculated by noting that the total distribution must be normalised, i.e.
i [P (i, +) + P (i, −)] = 1 and that J is imposed by the boundary conditions. For example, closed boundaries require that J = 0, while for periodic boundaries we have the additional constraint that P (L+1, ±) = P (1, ±). Note that the probability densities for left and right moving particles do not have to be the same, and, in general, will not be.
The probability of a given configuration P ( n + i , n − i ), is given by where N is the total number of particles in the system. If we call the average number of right or left going particles on a site ρ ± i = N P (i, ±), then in the limit where N → ∞ and L → ∞, so that P (i, ±) → 0, but ρ ± i remains finite, the probability of a configuration is given by 3.1.2. Continuous limit Since we ultimately want to compare the run-and-tumble on lattice with its off lattice counterpart, let us first take the continuum limit of the master equation. Explicitly introducing the lattice spacing a and defining x k = ka, the master equation (2) reads We are interested in cases where the typical run length is much longer that the lattice spacing so that d ± α ± . Furthermore, when a → 0 while v ± (x) = d ± (x)a remains finite, one gets, at leading order, which is exactly the master equation for run-and-tumble bacteria considered previously off lattice [12,13]. Following the same path as there would lead to a Langevin equation for the density of a large but finite number of bacteriȧ where with α = (α + + α − )/2, v = (v + + v − )/2. It was shown that (15) captures the steady state of the off lattice model exactly [13] and in section 4 we will show that it also describes the large scale behaviour of run-and-tumble bacteria on lattice. The condition for (15) to derive from an effective free energy is that there exists an excess free energy which can be solved to give as long as, for periodic boundary conditions, The total free energy is then given by Note that it is not always possible to write a free energy of this form for non-interacting particles in higher dimensions, nor, in general, for interacting systems. As we can see, there is no gradient term in this expression. This is due to the fact that when deriving the continuum limit, terms of order a and higher are neglected. If equation (15) leads to large gradients in the density, these higher order terms should be included and may alter the result; these terms control the surface tension, for example. These higher order terms could also violate the condition (17).
We will now illustrate the steady state predictions of equations (8) (9) and (19) for a few simple cases and show the difference between the results predicted by the lattice model and by the continuum approximation.

Examples
Where not stated otherwise the simulations we present here use reflecting boundary conditions; if a particle tries to jump off one end of the lattice it is, instead, kept where it is but turned around. All simulations are performed with continuous-time Monte Carlo algorithms.
Position Dependent Rates with Closed Boundaries First, we consider the case of a position dependent, but isotropic, jump rate and a constant tumbling rate. As a simple example we use a top-hat function for jump rate such that d ± i = 1+10 θ(i−150) θ(350−i), where θ(x) is the Heaviside step function. Both the continuum and lattice theory predict that the average occupancy should be inversely proportional to the velocity, ; The results of the simulations and both predictions are shown in figure 2 (main). In contrast to the jump rate, simply making the tumble rate depend on position but maintaining isotropy has no effect on the predicted distribution. Note that the free energy in equations (18) and (19) has no dependency on α for isotopic rates. Using Direction Dependent Rates with Closed Boundaries In many physical situations, however, bacteria do not move unbiasedly but are affected by their external conditions. This may be due, for example, to sedimentation due to gravity, where there is an asymmetry in jump rates between left and right (or up and down) moving particles. Another case of interest may be anisotropic tumble rates. Bacteria undergoing chemotaxis often vary their tumble rate dependent on whether they are travelling up or down a chemical gradient. Though a simple asymmetry in tumble rate does not fully capture this behaviour, we do see particles preferentially move in the direction of a lower tumble rate, as would be expected.
These two cases show qualitatively the same behaviour, with an exponential decay in the unfavoured direction. From equations (8) and (9), the probabilities for left or right going particles on lattice are while in the continuum case the probability density is given by For expositional simplicity we consider homogeneous rates here. To examine the difference between these two cases, consider, for example, the case of sedimentation, where d ± = d 0 (1 ± ) and α ± = α 0 . The decay constant in the continuum limit is then  Figure 3. Steady state probability distribution for constant tumble rates, α ± i = 1 and jump rates d ± i = 10 ∓ 1. Simulation data are shown in red, the lattice prediction in blue and the continuum prediction in pink. Data collected from 10 7 particles at t = 2000 on a lattice of 500 sites. and the lattice decay constant is given by We can see then that the two decay lengths will be equal if the jump rate is much larger than the tumble rate, i.e. for large average run lengths. Both decay constants tend to zero as the asymmetry disappears, → 0, but the ratio λ latt /λ cont remains finite. In our simulations we use d ± i = 10 ∓ 1 and α 0 = 1, so the drift velocity, the external bias, is much less than the self-propelled speed, i.e.
The continuum theory predicts the distribution to be ρ(x) ∝ exp(−x/99), while the lattice theory predicts P (i) ∝ (207/209) i . Both predictions are shown, along with the simulation data in figure 3. The ratio of decay constants is Note that in equations (23) and (24) one constant multiplies lattice position, i, and the other the continuum position, x, hence the factor of a. We see that the difference between the lattice and continuum results vanishes in the infinite run length limit, d 0 /α 0 → ∞, unless = 1, in which case both decay constants diverge.

Periodic Boundary Conditions
We can also calculate the expected probability distribution for periodic boundary conditions. In this case our calculation on lattice is slightly more complicated as we do not know the current a priori, but must determine it through the conditions P (L + 1, ±) = P (1, ±) and i [P (i, +) + P (i, −)] = 1. We can write P (i, where C 1 (i) and C 2 (i) can be read from equation (8). Then from the periodicity of the system we can write the current as The probability distributions at steady state for rates given by equations (28). The distribution for left moving particles is shown in blue and for right moving particles in red. The points show data from stochastic simulations and the solid lines the theoretical prediction. Data from 2, 000, 000 particles at t = 200 on a lattice of 200 sites.
We can then use the normalisation of the distribution to determine P (1, +) as As an example see figure 4 where we consider the case where the jump and tumble rates are We omit the exact forms for the probability distributions as these do not reduce to a compact form and are not enlightening in themselves.

A Zero-Range Interaction
Though there is no generic solution for the steady state of interacting run-and-tumble particles, there are limiting cases that can be solved exactly, most simply a zero range interaction. We define the number of left and right going particles on a site as, respectively, n − i and n + i . The total occupancy of a site is then n i = n + i + n − i . The simplest interaction we can add is a zero-range interaction, where the rate for a particle to jump from site i with occupancy (n and is a function of the number of particles at the departure site but not dependent on the number of particles at the arrival site. With this addition we can now see more complex behaviour and non-trivial steady states, even for homogeneous and isotropic jump and tumble rates, but can still, under certain conditions presented below, calculate the stationary probability distribution exactly. As is standard for zero-range processes (ZRPs) [18], we begin by assuming there exists a factorised form for the steady state probability distribution of the form This ansatz can then be substituted into the master equation Then, for periodic boundary conditions, one way we may choose to solve this equation is to separately balance the fluxes for right moving particles entering and exiting each site, left moving particles entering and exiting each site and particles tumbling between species on each site. We then arrive at three sufficient conditions on the allowed rates for such a factorised form to exist: in which c and c are arbitrary constants. The first two of these conditions are the same as Evans and Hanney found for their two species model without transmutation [18,20], while the third is due to balancing the tumbling. Note that, in principle, there may be other ways in which we can balance these terms which could arrive at different conditions on the rates. Putting the three conditions (31)-(33) together and eliminating the factors g i (n + , n − ) yields two constraints on our choice of rates: One natural, but again not necessary, way to fulfil these conditions is to take for the jump rates. That is, for the rate at which a left or right oriented particle moves to be a product of a function of the number of particles oriented in that direction, and a function of the total number of particles on a site. Both functions can vary from site to site; the first can also depend on the particle species, but the second must be the same for both. A sufficient conditions on the tumble rates is then to take The functions u + i (n + ) and u − i (n − ) are the same as in (36); A i (r) is a new, unconstrained, function that appears in both rates.
Up to a constant that can be subsumed into the normalisation, the marginals are given as where γ = c/c and n = n + + n − . We can then re-write the probability of a given configuration as where Z is a normalisation, and define an effective single site free energy ). Note that this is independent of α ± i (n + , n − ); the way in which we chose to solve the master equation does not lead to factorised steady states for which the distribution can depend on the tumbling rates. As we saw that asymmetric tumble rates could affect the equilibrium distribution for the non-interacting case, we might suppose there are other solutions for the zero-range process which admit distributions dependent on the tumble rates. Whether or not these allow for factorised steady states remains to be determined.
To foreshadow the finite range interaction we will examine in section 4, and to mimic the situation where an increase in density decreases the particles motility (as, for example, they get in each other's way) we now consider the following particular form of the steady state for this two-species ZRP for jump rates and tumble rate α ± i = α. That is, the tumble rate is constant per particle and the jump rate decreases linearly as density increases until reaching a constant rate of v 0 /n m at n + + n − = n m − 1. In this case the effective free energy is double welled and the system separates into isolated sites of high and low density. The relative numbers of high and low density sites to which the system first separates are initial condition dependent. The system then relaxes via a series of evaporations and condensations towards a fixed steady state. This is reminiscent of what happens for single-species zero-range processes with similar jump rates. Where the jump rates remain finite even for very large occupancies the normal condensation is arrested [21]. The free energy is shown in figure 5 along with a typical snapshot of the system.
To this behaviour we can then add a drift term to simulate sedimentation by biasing the jump rates in one direction and applying closed boundary conditions. We see all the high density sites collect at one end of the lattice and the low density sites at the other, see figure 6.

Fluctuating Hydrodynamics for Interacting Bacteria
We now turn to the more general case of N interacting bacteria with interactions that are not limited to being on-site. Specifically, we allow the jump and tumble rates to depend on the occupation numbers of each lattice site so that wheren ± i is a coarse-grained occupancy that depends linearly but non-locally on the occupancies of the whole lattice, In general the coarse graining kernel K ± i−j could also be a function of lattice position though here we do not consider that situation, where the manner in which the density is felt by the particles varies with position. We aim at describing large scale behaviour, i.e. on a colony size, and so to derive a fluctuating hydrodynamics. In addition, this will allow us to compare again the phenomenology on and off lattice and to look for cases in which there is a "free-energy" like description, and for which we can thus characterise the steady state. We follow a field theoretic approach to derive a continuum Langevin equation for the system, from which we can deduce the appropriate Fokker-Planck equation and the steady-state distribution.

Field Theory for Non-interacting Particles
To illustrate the technique we shall use to construct the fluctuating hydrodynamics for the full interacting case, let us first handle the non-interacting case and re-derive equation (15).
Field theoretic representations of lattice gases using bosonic coherent states were established in the 1970s following Doi and Peliti [22,23]. The case where each site is limited to a single particle can be handled in some cases in this formalism [24], while more general finite occupancies could be handled using spin coherent states [25]. Alternatively, probabilistic approaches from mathematical physics have also been used [26][27][28]. Here we use an alternative derivation, based on an approachà la Jansen and De Dominicis [29,30] transposed in the context of the master equation. This is very similar to the generating function approach used by Biroli and Lefevre [31]. Beginning with a process discrete in both time and space one writes the probability of a trajectory as where n ± i (t j ) is the number of right (+) of left (−) moving particles at site i at time t j , the J i (t j ) ± are the changes in the number of each type of particle at each site at each time step and the bold faced J denotes the average is over all J's. Re-writing the Dirac delta functions using imaginary Fourier representations this can be written as where it should be noted that the conjugate fieldsn ± i are imaginary. The average over the J's can then be calculated explicitly from the dynamics. Specifically, a right moving particle can jump from site i to site i+1 at time step j with probability n + i (t j )d + i dt, where dt is the duration of a time step. The corresponding values of the J's are J i (t j ) = −1 and J i+1 (t j ) = 1. Calculating all other moves and the probability that nothing happens, which corresponds to all J ± i = 0, we can write e −n + As this is of the form 1+kdt we can approximate it as exp(kdt) and write the probability for the trajectory as We can then take a continuous time limit and make the substitutions The probability of a trajectory can then be written where the action S is given by Note that generic changes of variables in (49) will result in Jacobians. If these do not depend on the fields, n ± i andn ± i , they can be subsumed into the normalisation of the path integral but they must be handled with care otherwise.
We further simplify by considering symmetric, constant, rates d + i = d − i = d and α + i = α − i = α; the more general case causes little conceptual difficulty but is considerably more cumbersome as an illustration. Let us then introduce the new variables The action can then be written as The continuum limit can be taken by explicitly introducing the lattice spacing a and making the substitutions where ∇ i is the discrete gradient, e.g. ∇ i ρ i = ρ i+1 − ρ i , and is the system length.
After Taylor expanding the action in powers of the lattice spacing and taking a diffusive rescaling of time and space, see Appendix A, we find that the fluctuating hydrodynamic action is given by which is invariant under a further diffusive rescaling of space and time.
Going back to the definition of the probability (49), we can then work backwards to recover a continuity equation for ρ from our action [25]. Starting from one can remove the quadratic termĴ 2 by introducing a new field η(x, t) via a Hubbard-Stratonovich transformation so that where the new action now reads (after some integration by parts) The integral overρ andĴ then leads to where the delta functions impose the two dynamic field equationṡ Given its weight in (58), η(x, t) is a Gaussian white noise: This is consistent with the calculation off-lattice for non-interacting, homogeneous and isotropic systems and validates the results obtained previously, see equation (15).

Fluctuating Hydrodynamics and Large Deviation Functions
Before going any further, let us make a brief detour to consider the connection with the standard fluctuating hydrodynamics approach considered in the mathematics literature [32]. Let us first note that from the definition of the continuum limit, one has The integral of the density field is thus an extensive variable. On the other hand, the density field considered by mathematicians is often defined by and satisfies the normalisation condition To make the connection between the two approaches, it is thus natural to rescale our density term to make the extensivity apparent:ρ → ρ. To ensure that the conservation equation still has the formρ = −∇J, one must also rescale the current field J → J. Before introducing the η(x, t) field, the action thus reads One can again introduce the noise field and integrate over the conjugates fieldsρ andĴ to get Interestingly, the fields are now all intensive, and the smallness of the noise does not come from a √ ρ versus ρ noise prefactor, but from its explicit variance, read in the Gaussian weight: This is the usual fluctuating hydrodynamics, as considered, for instance, in [32]. In the large size limit, the first order correction to the deterministic equation in a 1/ expansion is given by the addition of the noise term √ 2Dρ. This noise is typically of order 1/ √ , i.e. trajectories of probability of order 1 have η 2 of order 1. Large deviations correspond to trajectories where the noise can be of order one. They yield probabilities of order O(exp(− )), and are described by the fluctuating hydrodynamics constructed here.

Interacting Particles
Consider now the case of interacting particles where the jump and tumble rates depend on the occupation numbers of each lattice site. Our velocity is then modified to and the tumble rate to withρ ± (x) given by an integral over the densitȳ Following the same path as that followed in section 4.1 for the non-interacting case, one gets for the action where v = (v + + v − )/2 and α = (α + + α − )/2. The factor of in the final term implies that, for the diffusive scaling to hold, at a scale , the asymmetry must be of order 1/ . This is reminiscent of the ASEP, where if the bias is much smaller than 1/ √ the diffusive scaling holds (Edwards-Wilkinson universality class), as in the symmetric exclusion process, but for larger asymmetries the dynamic exponent z is the same as Kardar-Parisi-Zhang scaling [33]. Integrating overρ andĴ now yields the set of field equationṡ with This formalism can now be used to analyse the effect of interactions on the large scale behaviour of a system of run-and-tumble particles.

Crowding Interactions
Having set up and utilised a field theory apparatus to derive a fluctuating hydrodynamics for a general linear dependence on density in the jump rates, let us turn now to a specific class of interactions. In particular we consider a crowding interaction, where the velocity of the bacteria decreases with increasing density, in which case we expect to see our system separate into regions of high and low density as particles become trapped in regions of high density [13].
In general we expect to see qualitatively similar behaviour for any choice of v(ρ) which decrease sufficiently quickly towards a finite non-zero velocity at high densities. In the following we use as we did for the exactly solvable zero-range process.

Zero-Range Interactions
This approach can describe many types of interaction, in particular let us now consider the zero range interaction we met in section 3.2 in the context of the fluctuating hydrodynamics we developed in section 4. This provides us with a benchmark to check that our fluctuating hydrodynamics is consistent with the exact results we obtained previously.
For the zero-range process, where the velocity depends only on the occupation at the departure site, our kernel is given by K ± i = δ i0 and in continuumρ ± (x) = ρ(x). For simplicity, and to compare with our previous results we shall take α ± (ρ) = α. This simplifies equation (71) considerably and, indeed, guarantees that α + v − − α − v + = 0. From section 4.2 we know that the fluctuating hydrodynamics describing the run-andtumble bacteria with velocity given in equation (73) are given bẏ The corresponding Fokker-Planck equation is given bẏ Note that a term ∇ δ δρ [D(ρ)] could be present, but vanishes for symmetry reasons (See [13,34] whose solution is For ρ ≥ ρ m the free energy density f (ρ(x)) is given by which corresponds precisely to the free energy calculated exactly in section 3.2 for the total occupancy and in the continuum limit. An example of this free energy for one choice of parameters is shown in figure 5 (left).

Finite Range Interactions
We saw in section 5.1 that our fluctuating hydrodynamics admit a free energy and accurately reproduced the exact results for the zero-range process. Let us now extend that analysis to a system with finite range interactions. As before we now take the coarse grained densityρ ± to be given bȳ For smooth profiles, we hope that the differences betweenρ ± and ρ are small so we can treat the free energy in equation (77) as a mean field theory. This way we can still use the free energy we derived in section 5.1 to predict the coexistence densities and instability to spinodal decomposition. The finite range nature of the interactions will introduce correlations between sites which we hope will manifest only via surface tension [13]. We hope that this surface tension will only effect a clustering of the high and low density sites without affecting the coexistence densities. If the mean field theory captures the picture correctly, spinodal decomposition occurs whenever the second derivative of the free energy density is negative, i.e. for 1 ρ Thus, whenever ρ m > ρ > ρ m /2, the system should be unstable with respect to spinodal decomposition.

Isotropic Kernels
If we use an isotropic kernel to calculateρ ± in our jump rates we do indeed recover results consistent with the zero-range free energy. That is, if we take K + (x − y) = K − (x − y) our simulations match the free energy predictions. In particular, we have worked with a Gaussian kernel, where K ± (x) = exp(−x 2 /k)/Z, with Z a normalisation and k a parameter to control the range of the interaction. The results of simulations using this kernel are shown in figure 7 along with the predicted average high and low densities from the free energy. To calculate the coexisting densities of the high and low density sites we form a double tangent construction on the free energy [35]. As predicted the finite range nature of the interactions effectively creates a surface tension, but does not significantly alter the coexistence densities. Further, as expected we see no dependence in the steady state on our choice of v 0 and α.

Anisotropic Kernels
For anisotropic coarse graining kernels, however, the situation is more complex. One simple and natural way to introduce an anisotropic kernel is to account for the finite volume of bacteria by stating that there can be at most n m bacteria on each lattice site and taking the occupancy at the arrival site as our n in equation (43). This forms a generalisation of the partial exclusion process [25,36] and results in the jump rates In this case the effective free energy is limited to a region where ρ < ρ m and this section is sketched in figure 8. A double tangent construction amounts to looking for a density ρ low such that the tangent of the free energy density at this point meets with the free energy density at ρ = ρ m , as can be seen from inspection of figure 8. This amounts to finding ρ such that which can be solved numerically and yields for the low density ρ low that coexists with ρ m ρ low ρ m = .203188.
For a total average density larger that ρ low , the stable phase should thus be a combination of two phases, one with ρ = ρ low and one with ρ = ρ m , the ratio of the amounts of the two phases being set by constraint on the global mass. Interestingly, although the theory correctly predicts a change from a flat profile to phase separation, on examining the results of simulations of the underlying lattice system we found that the densities into which the system separates do not correspond to those predicted by the continuum theory. Indeed, while the continuum theory had no dependence on the tumble rate α or the coefficient of the jump rate v 0 , the simulations for an anisotropic kernel showed a strong dependence on the ratio of these two parameters. The lower and upper densities both varied considerably with the average run length r ≡ v 0 /α, as shown in figure 9, and below r = 4 we see no separation at all. This discrepancy is not limited to the particular choice of anisotropic kernel we use as illustration above and is general to any anisotropic choice of K ± (x). We have also conducted simulations with smooth anisotropic kernels without a hard limit on the  Figure 9. Snapshots of typical density profiles for an average run length of 100 sites (red), 20 sites (green) and 10 sites (blue). Data recorded at t = 1000 with α = 1, n m = 100 and from 10, 000 particles.
number of particles per site, and saw exactly the same qualitative effect. Thus we have seen that the relevant factor is indeed the isotropy of the kernel but although the origin of the difference has been established, a comprehensive explanation for the variation between isotropic and anisotropic kernels is yet to be formulated.
Note that it is not in itself surprising that a fluctuating hydrodynamics developed to describe smooth profiles fails to quantitatively amount for the coexisting densities of profiles alternating high and low densities. One of the reasons our fluctuating hydrodynamics work so well for the isotropic case is that we always consider large occupancies on each lattice site. This means that the model is close to mean-field in the same sense as the large spin limit of a spin chain is well described by a continuous spin chain [25]. Smaller occupancies would lead to quantitative differences between the predicted coexisting densities and those predicted from the fluctuating hydrodynamics, even for isotropic kernels. Furthermore, the Ito drift that was neglected in equation (75) for symmetry reasons would not vanish for the anisotropic case. In fact even for the offlattice model, the fluctuating hydrodynamics developed previously [13] is only valid for isotropic kernel and the quantitative mismatch between the fluctuating hydrodynamics and the simulations on lattice for anisotropic kernels are thus not that surprising. We nevertheless now try to shed some light on their origin.

Stability Analysis
One way we can analyse the difference between the isotropic and anisotropic interaction kernels is to examine the dynamic stability of the two systems. We consider a one dimensional system evolving from a homogeneous state under a small perturbation and determine whether the system is dynamically stable or unstable, whether the perturbations will, on average, grow or shrink.
One possibility is that the discrepancies we saw in figure 9 between the theory and simulations arise from the assumptions behind the diffusive limit taken in the field theory, see section 4 and Appendix A. We therefore start from the continuum microscopic mean field equations for homogeneous and isotropic rates, i.e. after the continuum limit has been taken but before the diffusive limit, We expand around a flat profile and Fourier transform. We take ρ ± (x) = ρ 0 /2 + q δ ± q exp(i q x) and arrive aṫ We can re-write these two equations in matrix form aṡ The eigenvalues of this matrix, which will then tell us whether the flat profile is stable or unstable to small perturbations, are It is clear that one of these eigenvalues will always be negative while the other is negative for ρ 0 < ρ m /2 and positive for ρ 0 > ρ m /2. Hence a homogeneously flat profile is stable when the average total density is less than half the maximum density and unstable above that, with no dependency on run length. This corresponds precisely with the stability predicted by the free energy derived in section 5.2. That stability analysis was derived from a free energy which considered only the total density and was itself calculated only after assuming a diffusive scaling. That the diffusive scaling does not alter the criterion for instability justifies taking that limit and implies that the discrepancy between our lattice simulations and continuum free energy arises from another factor. We turn, then, to consider the dynamic stability of the lattice dynamics directly. Beginning with the mean field equations for the anisotropic partial exclusion process, we expand around a flat profile, taking n ± k = n 0 + q δ ± q exp(i q k). After some algebra, detailed in Appendix B, we arrive at a condition for there to exist positive eigenvalues, i.e. for a flat profile to be unstable to small perturbations. Specifically we see instability whenever the run length r = d/α satisfies the following inequality Run Length I II Figure 10. A flat profile is stable when in region I, but unstable in region II, for a system with exclusion and homogeneous, isotropic jump and tumble rates. The x-axis is the fractional density, i.e. n/n m .
Graphing this we can see that a flat profile will be stable in region I of figure 10 and unstable in region II. We see that for short run lengths the range of densities in which the system will spinodally decompose is restricted and at run lengths below 4 sites there is no separation at all. Conversely, in the limit that the run length tends to infinity, i.e. where we effectively have two totally asymmetric partial exclusion processes on the same lattice, the system is unstable for any density between n 0 nm = 0.5 and n 0 nm = 1. This instability is in accordance with that seen in our simulations.
In our simulations we found that when we replaced the anisotropic kernel in the interaction terms with an isotropic one we recovered the density profiles predicted by the continuous theory. We can also analyse the effect of an isotropic density kernel on the dynamic stability.
Consider now the microscopic mean field equationṡ In general K ± j can take any values, we enforce only that they are both normalised, i.e. that j K ± j = 1. Relaxing this constraint would effectively re-normalise the maximum density. This more general interaction reduces to the simple exclusion case if we take K ± j = δ ±1,j and to the zero-range case if we take K ± j = δ 0,j . For an isotropic kernel, where K ± i = K i , when we expand around a flat profile we find that there exists a q such that λ + (q) is greater than 0 if and only if n m > n 0 > 0.5 n m , see Appendix B for details, which matches the condition we derived from our continuum free energy. Thus for an isotropic interaction kernel we recover the continuum stability result while an anisotropic kernel will, in general, not produce the same result. It thus seems the error comes from the continuum limit itself, which is not valid for anisotropic kernels. While the stability analysis accounts qualitatively for the difference between isotropic and anisotropic kernels a theory which quantitatively accounts for the differences in coexistence densities at steady state remains to be constructed. Figure 11. Snapshots of the two dimensional system using the partial asymmetric exclusion and parameters n m = 20, d 0 = 10 and α = 1. Sites with a density greater than 0.6n m are marked in red, while those with a lower density are left blank. Left: allowing only nearest neighbour hops and no diagonal movement. Right: allowing for diagonal hops. Both simulations performed using 400, 000 particles and recorded at t = 2500.

Finite Range Interactions in 2D
Most cases of physical interest require a model in more than one dimension; it is therefore natural to extend our analysis to higher dimensions. This can be done on lattice relatively easily. On a square lattice in two dimensions we allow the particles to jump between nearest neighbours. We therefore consider 4 types of particle now instead of 2. We find that in two dimensions, the behaviour of the run-and-tumble crowding model is qualitatively the same as in one. The system separates into regions of high and low density, where those co-existent densities are given by the same free energy as in one dimension. The field theoretic approach developed in section 4 indeed generalizes straightforwardly to higher dimensions and yields the same fluctuating hydrodynamics. We find, however, that allowing only nearest neighbour hopping results in an unrealistic surface tension because of the anistropy of the lattice [37,38]; the regions of high and low density form elongated, and thus anisotropic, domains, see figure 11 (left).
To correct this unphysical characteristic of our model we extend our dynamics to allow next to nearest hopping along diagonal directions and consider 8 species of particles, where the jump rates in the diagonal directions are scaled by a factor of √ 2. The droplets then coarsen into more realistic curved domains, see figure 11 (right).
The stability conditions remain the same as do the coexistence densities and the discrepancy between isotropic and anisotropic kernels. We examined simulations with both the anisotropic partial exclusion kernel, and the isotropic, Gaussian kernel on lattices of 500x500 sites for densities above the spinodal point. The systems were seen to separate into droplets of higher and lower density which then coarsened into discrete, contiguous domains, see figure 11.
We measured the coarsening of these domains and found them to scale as t 1/3 , as we would expect for conserved model B type dynamics [39]. An approximate measurement of the size of the domains was calculated by randomly sampling the system and measuring the horizontal and vertical size of the encountered droplet at that point. Mathematically, we define where n * is an arbitrary cutoff to distinguish the two domains but ignore random fluctuations. Computations were run with a number of choices for n * and the particular choice of cutoff was found to have no significant effect on the results. We calculate the vertical size in an analogous fashion and average the lengths over a large number of points on the lattice. Though this does not give an exact measure of the droplet size it is sufficient to show the scaling of the domain size with time whilst being quick to calculate numerically. Using this procedure we determine that the domains increase in size with an exponent of approximately 1/3, i.e. L x (t ) = (t /t) 1/3 L x (t). The two-point, connected, equal time correlation function C(j, t) = n i (t) n i+j (t) − n i (t) n i+j (t) was also calculated numerically from the data and fit reasonably with a re-scaling C(x, t) = C(x/a 1/3 , t a). The data for both these measurements can be seen in figure 12.
As can be seen from the figures the choice of kernel does not change the coarsening exponent, only the relative speed of coarsening, with interactions over a larger number of sites taking longer to reach a steady state than those with shorter ranges.

Conclusion
In this work we have presented several classes of lattice models based on the run-andtumble dynamics of certain species of bacteria, notably Escherichia coli. We calculated the exact steady state probability distributions for both inhomogeneous, anisotropic, non-interacting and zero-range interaction models. For more general types of interaction we used a field theoretic approach to derive the continuum fluctuating hydrodynamics and from there derived a mapping to a free energy like functional describing the steady state profile for a crowding interaction. We analysed the linear stability of both the continuum and lattice microscopic mean field equations and isolated a condition on the coarse graining we employed in our interaction terms for the zero-range free energy to work as a mean field theory.
Our work builds on earlier treatments of run-and-tumble bacteria where interactions between bacteria were not included [12,40]. It provides a lattice counterpart to prior continuum approaches and qualitatively reproduces results obtained off lattice, where a similar, though not exactly equivalent, density dependence was considered [13]. Our approach on lattice provides a microscopic justification for manner in which this density dependence manifests; previously it was added in an ad-hoc manner after the diffusive approximation had been taken where in this work the dependence is intrinsic from the microscopic definition of the dynamics. That we produce qualitatively similar results justifies the way in which the dependence was inserted there. Our work also reveals a condition on the way in which the coarse graining in this density dependence must be performed for free energy mapping to work; it must be taken isotropically, so that particles moving in any direction feel the same effective local density.
This work provides a means to simulate microscopic run-and-tumble dynamics efficiently, which is particularly important in two or more dimensions where microscopic simulations off lattice are very computationally expensive. It illustrates potential hazards in comparing lattice simulations and continuum theoretical predictions and offers some insight into how to avoid those problems by carefully choosing how to implement non-local interactions in the lattice dynamics. Note that early studies of off-lattice run-and-tumble dynamics actually had to discretise space to do simulation for position dependent swimming speed in order to compare with their theoretical predictions [40]. We hope the techniques we presented in this paper will develop the use of lattice simulations in the bacterial context much further, for example, to study pattern formation.
We hope that the lattice approach will provide new exact results for run-and-tumble dynamics in higher dimensions, generalizing previous studies to more general external potentials for sedimentation and trapping [14].
Last, there are some open question regarding the coexistence densities for interacting bacteria. Although our fluctuating hydrodynamics correctly predicts the existence of phase separation, the coexistence densities are accurately predicted only for isotropic kernels and large lattice occupancies. Quantitative predictions beyond this case are yet to be derived. More generally, explicit two body interactions and more general non-linear interaction kernels have not been looked at yet. These are a few of the outstanding questions among the many concerning this new and interesting class of models.
If we explicitly takeĴ to scale asĴ →Ĵ −β and J → J −δ we get Now, we need the coefficient of each term to be of order 1 or smaller so that no terms blow up so However, as we do not want to simply be left withρ = 0 we need α − δ = 0 and as we also want to retain a noise, which corresponds to theĴ 2 term, we also require that α − 2β = 0. That leaves our action Which tells us that β ≤ 1, β ≥ 1/3 and β ≥ 1 which imply that β = 1 and hence α = 2 and δ = 2. Injecting these scalings back into the action gives where the macroscopic observation time τ = T / 2 is supposed to be of order 1. The term inĴJ is thus irrelevant and we can check that the hydrodynamic action is invariant under further diffusive scaling. Note that the scaling of the fields considered here is arbitrary and we could choose to look at currents J,Ĵ larger than 1/ 2 , 1/ . This would correspond to trajectories whose probabilities are smaller than exp(− ), which are not correctly described by fluctuating hydrodynamics and large deviations. One can also check that under this rescaling, the action S 1 stays of order 1 and aS 1 is thus, indeed, negligible. For one isolated bacterium the run-and-tumble dynamics is a variant of a random walk and is known to be diffusive at large scales [12,13]. It is therefore not surprising that we find α = 2. In the presence of interactions (as in section 4.2), a uniform density profile of bacteria will continue to exhibit diffusive behaviour; the interactions will simply rescale the diffusivity. If interactions cause the profile to become unstable, however, the model can nevertheless give rise to length scales which grow in a non-diffusive manner, as for example for coarsening (see section 5.3).