Schr\"odinger's ants: A continuous description of Kirman's recruitment model

We show how the approach to equilibrium in Kirman's ants model can be fully characterized in terms of the spectrum of a Schr\"odinger equation with a P\"oschl-Teller ($\tan^2$) potential. Among other interesting properties, we have found that in the bimodal phase where ants visit mostly one food site at a time, the switch time between the two sources only depends on the ``spontaneous conversion"rate and not on the recruitment rate. More complicated correlation functions can be computed exactly, and involve higher and higher eigenvalues and eigenfunctions of the Schr\"odinger operator, which can be expressed in terms of hypergeometric functions.


INTRODUCTION
Kirman's ant model [1] undoubtedly stands among some of the most inspiring toy models in the behavioral economics literature. While initially inspired by the experiment described below, its conclusions have implications much beyond collective animal behaviour, and has been used to model shifts in sentiment of economic agents, trend reversal in financial markets, herding and social influence, etc. Kirman's model is also akin to another famous model in population dynamics with competing species: the Moran model [2].
Several decades ago entomologists were puzzled by the following observation [3,4]. Ants, faced with two identical and inexhaustible food sources A and B, tend to concentrate on one of them for a while, but occasionally switch to the other. Such intermittent herding behavior is also observed in humans choosing between equivalent restaurants [5], or in financial markets [6][7][8] consistent with large endogenous fluctuations. Clearly the asymmetric exploitation observed in ants does not seem to correspond to the equilibrium state of an isolated representative ant with rational expectations. The phenomenon is rather to be explained in terms of interactions between individual ants, or, as put by biologists, recruitment dynamics. To account for such intricate behavior, Kirman proposed a simple and insightful model [1] based on tandem recruitment that we now recall.
Consider N ants and denote by k(t) ∈ [0, N ] the number of ants feeding on source A at time t. When two ants meet, one of them converts the other with probability µ/N , but each ant may in addition change its own mind spontaneously with probability . Within such a simple setting, Kirman was able to show that, in the large N limit, the stationary state depends only on a parameter α := /µ. When α > 1 the distribution is unimodal, with a maximum at k = N /2, whereas for α < 1 the stationary distribution of k is bimodal, with maximum probability for k = 0 and k = N (corresponding to the situation observed in the experiments). It is remarkable that the interesting α < 1 regime can be obtained even for weakly persuasive agents (µ small) provided self-conversion is itself low enough.
The most important point is that in the α < 1 regime no one of the k states is, in itself, an equilibrium. Although the system can spend a long time at k = 0, N (local stationarity) these states cannot be considered as such: all the states are always revisited and there is no convergence to any particular state, discarding also the notion of multiple equilibria. Rather, there is perpetual change, and the system's natural endogenous dynamics is only in a statistical equilibrium. Most economic models focus on finding the equilibrium to which the system will finally converge, and the system may only be knocked off its path by large exogenous shocks.
Yet financial markets, and even larger economic and social systems, display a number of regular large switches (correlations, mood of investors etc.) which do not seem to be always driven by exogenous shocks. In Kirman's stylised setting such switches can be understood endogenously. Several extensions of the model have been proposed [8][9][10]. In particular, the original version of Kirman's model does take into account the heterogeneity in encounter probabilities induced by the topology of the social network; but one can easily (at least numerically!) modulate the probability of encounters according to their distance along such a network, for example restricting recruitments to nearest neighbours only.
In the present paper we present a continuous description of Kirman's ant model which notably allows us to derive the typical switching time, using classical methods from statistical physics and quantum mechanics.

I. MASTER EQUATION
As mentioned above the original model describes N ants faced with two identical food sources, with the relevant dynamical variable being k, the number of ants feeding on -say -source A. Each time step allows an ant to either switch randomly to the other food source with probability proportional to , or get recruited by another ant from the other food source with probability proportional to µ.
Defining the unit of time as the time required for all the ants to make a decision, leads to dt = 1/N as the infinitesimal time step. It is also clear that, to remain intensive in the large N limit, the probability to interact with another ant should be proportional to 1/N . Altogether, we may write a Master equation for the evolution of the probability P(k, t) that there are k ants feeding at source A at time t: where the transition rates are given by: Note that this specification only differs from Kirman's original one in the rescaling of the recruitment rate by N . With the notations of [1], 1 − δ = µ/N .

II. CONTINUOUS DESCRIPTION AND FOKKER-PLANCK EQUATION
Here we follow Kirman's original paper [1] and derive a proper continuous-time Fokker-Planck equation in the limit N → ∞.
We define the variable x = k N ∈ [0; 1] together with its probability density function f (x, t). Taking the continuous limit N → ∞ of Eq. (1) leads to the following Fokker-Planck equation [11]: the probability flux, see Appendix A for the details of the calculations and the first 1/N corrections. The conservation of the number of ants in the model is ensured by the condition J f (x, t) = 0 at the boundaries x = 0 and x = 1 at all times. Equation (3) corresponds to the following stochastic process for x: with η a Gaussian white noise with unit variance. One can note that while the drift term (1 − 2x) is maximal at the boundaries and tends to pull x towards 1/2, the noise term has the opposite effect. The diffusion constant is proportional to 2µx(1 − x) and is maximal at x = 1/2 and so tends to push the system away from x = 1/2. Note that this stochastic process is very similar to the Moran model of genetic population dynamics [2] -with the same diffusion term ∝ x(1 − x) -where the analogue of the number of ants at each food source is the proportion of genes from two competing alleles (A or B) [12]. The term corresponds to spontaneous mutations. When = 0, there is a non zero probability that the whole population becomes of type A or B after a finite time, corresponding to δ(x) or δ(1− x) contributions to f (x, t) with a time dependent weight, see [13], and [14] for a recent thorough discussion.
When > 0, one can check that the normalised stationary distribution f 0 (x), obtained by setting J f (x, t) = 0, writes: This result is the same as that obtained by Föllmer and Kirman in [1].
Upon looking at the behaviour of the solution, shown in Fig. 1, one can see that there is a clear transition in the behaviour of the model at α c = 1. For α > α c , the stationary density in Eq. (5) is maximal at x = 1/2, and the dynamics shows that x(t) fluctuates around 1/2, corresponding to a situation where the ants are, on average, evenly distributed across both food sources. For α < α c the density f 0 diverges at the boundaries. The top left panel in Fig. 1 shows that this corresponds to a very different picture, in which nearly all of the ants choose either one of the sources for a certain amount of time, until a noise-induced "avalanche" causes a switch over to the other source. It is also easy to check that in the absence of noise (and α → 0) the long-time stationary density is given by ], a situation discussed at length in [14].
Having this in mind, a natural question to ask is: Given a certain initial condition f (x, 0) = δ(x − x 0 ), how long does it take for the system to converge to the stationary state, or equivalently, how long does it take for the ants to switch from one source to the other in the α < 1 regime?

III. SCHRÖDINGER'S EQUATION AND GENERAL SOLUTION
Here we obtain a full dynamical solution in terms of the eigenvalues and eigenfunctions of a certain quantum mechanical Hamiltonian.
Using the Itô rule [15], one can see that introducing a change of variables ϕ(x) in Eq. (4) yields a noise term proportional to x(1 − x)ϕ (x), and so motivates a choice satisfying ϕ (x) = 1/ x(1 − x). We therefore define a new, more convenient, variable ϕ ∈ [−π/2, π/2] as: The corresponding Fokker-Planck equation for its probability density g(ϕ, t) writes: where the probability flux must now verify J g (±π/2, t) = 0 at all times. Setting again J g = 0 everywhere, one finds the normalized stationary solution: The advantage of this formulation in ϕ is that, in contrast with the former, the second order derivative term ∂ ϕϕ in Eq. (7) only depends on ϕ through g(ϕ, t). Standard techniques for the resolution of Fokker-Plank equations, see e.g. [11], motivate the introduction of a function Ψ such that: and Ψ(ϕ, t) → g 0 (ϕ) when t → ∞.
Combining Eqs. (7) and (9) one obtains a Schrödinger-like equation of the form [16]: where the Hamiltonian H is defined as: and with boundary conditions given by: We have left the µ parameter out of the Hamiltonian H in order to ease the comparison to the canonical form presented in [17,18]. The tan 2 term in Eq. (10) is known as the Pöschl-Teller potential [17], which was fully solved in the case β > 0 with boundary conditions Ψ(±π/2, t) = 0 in [18]. To be applicable to our framework, we shall verify that their solutions also satisfy Eq. (12) in the general case β > −1/2. The Hamiltonian H is Hermitian (contrarily to the Fokker-Planck operator) and has a discrete set of orthogonal eigenfunctions and eigenvalues, given by: where, splitting into even (n = 2k) and odd (n = 2k + 1) states: with 2 F 1 the ordinary hypergeometric function. 1 The coefficients A n are set such as to ensure normalisation, [−π/2,π/2] Ψ n Ψ m = δ n,m , and can be expressed as integrals of hypergeometric functions. Note that the parity of n also defines the parity of the function Ψ n with respect to the y-axis. One can then easily check that for all n (both even and odd): which, since β > −1/2, ensure that the boundary conditions given by Eq. (12) are satisfied. Noting that Ψ 0 = g 0 , the general solution of Eq. (10) then reads: with λ n given by the projections of the initial conditions on each mode n, namely λ n = . Further using Eq. (9), it is easy to see that the initial condition in turn translates into Ψ(ϕ, 0) = δ (ϕ − ϕ 0 )/ g 0 (ϕ). The full solution for g(ϕ, t) follows: with the orthogonality between Ψ 0 = g 0 and Ψ n for n > 1 ensuring that π/2 −π/2 dϕ g(ϕ, t) = π/2 −π/2 dϕ g 0 (ϕ) = 1, or equivalently for f (x, t): with: (see Appendix D 2 for an explicit expression). Equation (18) is the central result of the present communication.

IV. RELAXATION TOWARDS THE STATIONARY STATE
With the full dynamical solution of Eq. (18) at hand, one can see how long a system initially prepared at an initial value x 0 ≈ 0, for example, takes to explore the whole space. In other words, one can ask how much time τ is required to reach, say, x(τ) ≈ 1 with a reasonable probability. 1 Here, the function 2 F 1 takes the form of a polynomial: 2 u for any integer k. Since the stationary distribution f 0 has weight on the whole interval [0; 1], this time τ is none other than the relaxation time (or ergodic time) τ R required to converge to stationarity. Owing to the form of Eq. (18) this convergence is asymptotically exponential, with the slowest mode given by n = 1. Hence, we find: Perhaps surprisingly, this relaxation time depends only on the noise intensity , but not on the recruitment intensity µ. Since n = 1 corresponds to the slowest mode of the system, it also governs the collective "switch time" between the two food sources, A and B -see Fig. 2.
We have checked our prediction for the switching time numerically by running trajectories starting at x 0 = ∆x 1 and computing the probability (x(t) > 1 − ∆x). This quantity should converge to [1−∆x;1] f 0 at an exponential rate ∝ e −µ 1 t , which is in perfect agreement with our simulations, see Figure 3. Similarly, given an initial condition x 0 = 1/2 where the ants are initially distributed evenly between the two sources, one may ask how long it takes for all the ants to "decide" on concentrating on one of them. Since this condition is equivalent to ϕ 0 = 0, and since Ψ 1 is an odd function of ϕ, it follows that Ψ 1 (ϕ 0 ) = 0 in this case. The convergence to the stationary distribution is then controlled by the second mode, with a much shorter relaxation time given by: Directly applying tools from stochastic calculus on Eq. (4), one can obtain the following correlation functions (see Appendix E): where σ n (x) are polynomials of degree n that allow one to "diagonalize" the evolution of the correlations: See Appendix E for further details and Figure 3 for a comparison with numerical results.
where the exact expression of B n,m is given in Appendix D 3, Eqs. (D20) and (D21). Mind that B 0,m is the stationary value of moment [x m (t)] for all moments.

V. CONCLUSION
In this work, we have shown how that the approach to equilibrium in Kirman's ants model can be fully characterized in terms of the spectrum of relaxation times, itself computable as the eigenvalues of a Schrödinger equation with a Pöschl-Teller (tan 2 ) potential. Note that similar techniques have been recently applied to discuss the dynamics of wealth inequality in Ref. [19]. Among other interesting properties, we have found that in the bimodal phase where ants visit mostly one food site at a time, the switch time between the two sources only depends on the "spontaneous conversion" rate and not on the recruitment rate µ. This means that a single ant deciding on its own to explore an alternative food source can trigger an "avalanche" where the whole colony follows suit. More complicated correlation functions can be computed exactly, and involve higher and higher eigenvalues and eigenfunctions of the Schrödinger operator.
The possibility to solve exactly the dynamics of Kirman's model is of course intellectually satisfying. It is also important in view of the number of possible applications of such a model, recalled in the introduction, and which has reappeared recently in the context of self-fulfilling prophecies in a simple economic model [20] and in the empirical study of the dynamics of fishers seeking to exploit fishing zones with finite resources [21]. Our analytical approach can also be easily generalized to other models of genetic population dynamics, such as the general setting discussed in [14], as the change of variable we introduce always leads to a Schrödinger equation with a trigonometric potential provided the drift is linear in x. These equations may then be solved using known analytical tools [22].

We warmly thank Roger Farmer, Alan Kirman and Joachim Krug for fruitful discussions. This research was conducted within the Econophysics & Complex Systems Research Chair under the aegis of the Fondation du
Risque, a joint initiative by the Fondation de l'École polytechnique, l'École polytechnique and Capital Fund Management.

Appendix A: Derivation of the Fokker-Planck equation and stationary solution
We define the continuous distribution f (x, t) as: which amounts to replacing k N by x in Eqs. (1) and (2). In this case, and to leading order in 1 N , the term e.g. W (k + 1 → k)P(k + 1, t) reads: We proceed similarly for all terms in the right-hand side of Eq. (1), and Taylor-expand the left-hand side to leading order in the time variable, to obtain: where ∆ = 1 N for simplicity. We next Taylor-expand the right-hand side terms, such as e.g.
, to order ∆ for the terms with prefactor /∆ and to order ∆ 2 for the terms with prefactor µ/∆ 2 . Gathering everything, we obtain the Fokker-Planck equation: the same as given in Eq. (3). This equation can be written as , where J f is the probability flux, a function such that J f (x)∆ corresponds to the probability mass flowing from x + ∆ to x. To ensure the conservation of probability in [0; 1], we impose J f = 0 at the boundaries, meaning that no probability mass comes in or goes out during the dynamic evolution of the process. In other words, writing I f (t) = Recalling now that a Fokker-Planck equation of the form corresponds to the Itô stochastic differential equatioṅ where η is a brownian white noise, one readily recovers Eq. (4). Physically, the 0-flux boundary condition corresponds to a reflecting boundary condition: a "wall" that prevents x from getting out of [0; 1].

Determining the stationary solution
Looking for a stationary solution, one sets the right-hand side of Eq. (A4) to 0, looking to solve which, after direct integration, yields f 0 (x) ∝ (x(1 − x)) α−1 . Integrating for x ∈ [0; 1] allows one to find the normalisation constant in terms of the Beta function, or equivalently as a ratio of Gamma functions, to get Eq. (5).

Appendix B: Change of variables under an SDE
Obtaining Eq. (7) and understanding the rationale behind the change of variables of Eq. (6) is easier by starting from Eq. (4).
Imposing a change of variables x → ϕ(x) leads to a new stochastic differential equation for ϕ, which after applying the Itô rule for differentiation reads which is still difficult to interpret because of the dependence on x of the term in front of the white noise η.
Keeping instead the term of order ∆ given in Eq. (A5) leads first to the Langevin equatioṅ which leads to the change of variables where now |φ| ≤ arctan 1/ 2α∆ ≈ π 2 − 2 α∆, and naturally one can check that the definition of φ corresponds to ϕ as ∆ → 0, with φ ≈ ϕ − 2α∆ tan(ϕ) to leading order in ∆. The analysis in the limit N → ∞ therefore holds only in the limit tan(ϕ) N 2α . This new variable actually verifies the very same SDE, Eq. (B1), but with a different boundary.

Appendix C: Schrödinger from Fokker-Planck
The following is a common "trick" to transform a non-hermitian dynamic evolution coming from a Fokker-Planck equation with drift into a hermitian evolution determined by a Schrödinger equation. We start from a generic Fokker-Planck equation such as the one defined in Eq. (A6), but with constant b( y, t) = 1 and timeindependent drift, which we represent with the derivative of some function A, a( y, t) = −A ( y). The resulting Fokker-Planck equation reads and has a stationary solution that can be written as a Boltzmann distribution p 0 ( y) = e −A( y) /Z, where Z is a constant ensuring normalisation. We next introduce a function Ψ verifying p( y, t) = e −A( y)/2 / ZΨ( y, t). We can compute derivatives to find Adding these terms and simplifing, we find the following Schrödinger's equation for Ψ: where the Hamiltonian is here defined as Equation (10) simply uses this substituion, with dϕ tan ϕ = log cos ϕ playing the role of A( y) (up to a multiplicative constant).
which, written as such, leads us to introduce the function Applying the generalized Leibniz rule to compute the k − p − 1-th derivative of this function, we obtain directly that S 2k,2p = P (k−p−1) (1) = 0 in this case. A similar calculation can be done for S 2k+1,2p+1 , and it follows therefore that while on the other hand successive integration by parts gives Finally, gathering everything we get while replacing β → β + 1 gives the similar expression The final result follows, with I 2k,2p+1 = 0 allowing then for explicit computation of the dynamics of [x m (t)].

Appendix E: Stochastic calculus techniques
In this Appendix, we shall directly integrate stochastic differential equations describing the model to obtain information on the covariances of moments x n (t). We begin by looking at the covariance Cov(x(t + T ), x(T )).
This method can be extended to computing C n,k (t + T, T ) = Cov x(t + T ) n , x(T ) k . Applying Itô calculus as before, one can show that these functions satisfy the following ODE system: d ds C n,k (T + s, T ) = −µ n C n,k (T + s, T ) + µn(n − 1 + α)C n−2,k (T + s, T ).
Owing to its triangular structure, it can be diagonalized iteratively to find functions σ n , such that σ n (x) is a polynomial of degree n and that the covariances C σ n (T + s, T ) = Cov [σ n (x(T + s)), σ n (x(T ))] satisfy d ds C σ n (T + s, T ) = −µ n C σ n (T + s, T ).
These results can also be obtained directly from the eigenvalues and eigenfunctions of the Schrödinger problem.