of the number of objects on each site, which can eventually break the unitarity constraints. Nevertheless, we do not a priori specify a saturation mechanism at this level: Several options of implementing unitarity of the scattering amplitudes will be examined in the next section, and saturation, in the sense that the number of particles is prevented from exhibiting an exponential growth forever, will only be one of them.

Note that the choice of discrete transverse space and especially discrete $y$ is not very natural for a model that is meant to mimic QCD, in particular since discretizing $y$ obviously breaks boost invariance (which makes sure that the rapidity evolution can be shared arbitrarily between the two incoming hadrons without affecting the observables). However, this choice is dictated by technical reasons: We will need to be able to generate millions of events within some reasonable computer time.

Let us define the generating function for the probability of the different configurations as $$Z(y,\{x_k\})=\sum_{\{n_k\}} \left(\prod_k x_k^{n_k}\right) P(y,\{n_k\}),$$ where $P(y,\{n_k\})$ is the probability of having $\cdots n_1,\cdots,n_k,\cdots$ particles on sites $\cdots 1,\cdots,k,\cdots$ at rapidity $y$. This function contains all information about the statistics of the particle numbers on all sites. For example, the correlator 