A Mexican hat dance: clustering in Ricker-potential particle systems

Free of confinement (i.e., with $w\to\infty$) and in the absence of any degenerate initial positions, particles will settle into a state of uniform spacing, with each particle residing at the preferred distance $s$ from its neighbors. However, when confined, they exhibit highly nonuniform and rich behavior.

When confinement is present, we observe spontaneous organization of the particles depending on the choice of confinement parameter $w$ (note that confinement is stronger as $w$ decreases toward zero). Particles form “stacks”—states where multiple particles occupy the same spatial position—since their repulsion weakens as they get nearer to each other, but the number of particles in each stack can vary, and indeed the stacks exchange particles as $w$ changes. We demonstrate this spontaneous organization in Fig. 3, which shows the results of numerical investigation of the system’s stable equilibria.

Note that there is some slight asymmetry in Fig. 3 due to high dimensional multistability with various stack sizes; if we enforce symmetry, we see a picture of “core” behavior as shown in Fig. 4.

Also, note that these figures actually display $2$-dimensional projections of an $n$-dimensional bifurcation diagram, with each particle’s position occupying one dimension. However, the interchangeability symmetry of the particles allows the display of the whole population to serve as analogue for any single particle’s possible positions.

The stack-size information (encoded by color) can be further emphasized with a third dimension, as shown in the right panel of Fig. 4.

Thus we see even more clearly the pattern of stable behavior. For very small $w$ (strong confinement, i.e., a narrow parabolic well) all the particles stack up at the origin, but the population splits apart into two symmetric stacks¹¹1There is some multistability with regards to slight asymmetry of these two stacks; for example, sometimes the stacks are 63 and 65 particles, and the system is stable, although in that case their positions are not perfectly symmetric—the larger stack settles slightly closer to the origin. when $w$ passes a critical threshold we label $w_{0}$—we derive a formula for this value (Eq. (12)) in the “Large $n$ analysis” section. As $w$ continues to grow, the stacks drift apart and the origin becomes stable again, and we see particles “fall” inwards to settle there. Initially only a few particles stably rest there, but as the stacks continue to separate the central stack grows, until it becomes large enough to split into two in a manner that appears similar to its initial bifurcation at $w=w_{0}$.

There are many other, more complex equilibria possible, but for large $n$ the equilibrium diagram appears to become increasingly well characterized by the aforementioned behavior, as shown in Fig. 5. All the “major” bifurcations (birth of central stacks, splitting of central stacks) appear to happen at the same multiples of the critical parameter value $w_{0}$, thus with proper scaling of the $w$ axis (to match $w_{0}$) the diagrams appear increasingly similar to one another.

We will start our exploration with smaller, more tractable examples and then progress from there to the more general, large-$n$ cases.

We can solve for the equilibrium condition on the right particle by assuming the particles are at $\pm x^{*}$ (plugging in $2x^{*}$ as the distance in equation (3) and adding the background potential at position $x^{*}$):

At equilibrium this slope is zero, thus $x^{*}=0$ (corresponding to both particles stacked at the origin) or

So the bifurcation diagram (in $w,k,$ or $s$) is given by the implicit equation (5). This corresponds to the exact solution

where $W$ is the Lambert W function, defined as the solution to

Fig. 6 shows this solution overlaid on the empirical equilibrium diagram for two particles.

An approximation for $|x|\ll 0$ yields

where we have defined $w_{0,2}=\frac{s}{k}\sqrt{\frac{k-1}{2}}$ as the lowest critical value of $w$ with $n=2$. We derive the general-$n$ version (which we simply refer to as $w_{0}$ going forward) in the “Large-$n$ analysis” section (Eq. (12)).

Three- and four-particle systems have unique stable equilibria for all $w$, which are shown in Fig. 7; unfortunately these resist such easy exact-solution form as the $n=2$ case. Other equilibria exist for these systems (such as $1-2$ states²²2This notation indicates the partition state of the particles, i.e., the number of particles in each stack, in order. In this case, there are two visibly-distinguishable $1-2$ states since the two-particle stack can be above or below the single-particle stack, though each of those represents six equilibria in full state space due to permutations of the particles. in the $n=3$ case, $1-2-1$ and $1-3$ states in $n=4$, and the fully-stacked origin state at $w>w_{0}$), but none of them are stable. In Fig. 8 we show all such equilibrium positions using the MatCont analytical-continuation software package for Matlab (MatCont v7.3, see [dhooge2008new]) to track those unstable artificially-partitioned states as well as the stable state we actually see in regular numerical simulations.

The nonzero stable branches for $n=3$ (and $k=2,s=1$) are solutions to the implicit equation

which is the result of assuming symmetry (one particle at $0$ and the other two at $x^{*}$) and solving for nonzero $x^{*}$ solutions to the equilibrium condition on the particle at $x^{*}$. For large $w\gg w_{0}$ this relationship converges to $w^{2}=e^{2x^{2}}(4-4x^{2})^{-1}$ (or explicitly in $x$: $x^{2}=1-\frac{1}{2}W(\frac{1}{2}e^{2}/w^{2})$); for $w$ near $w_{0}$ ($w_{0}=\sqrt{3}/6$) the relation is well approximated by $x^{2}\approx 1/9-1/108\ w^{-2}$ instead.

In the 4-particle case, ($n=4,k=2,s=1$) we can find the implicit equation for the $2-2$ (i.e. two stacks of two particles each) state, namely

which is stable from $w=w_{0}$ ($=1/4$ in this case) to $w\approx 2.1w_{0}$ as seen in Fig. 7. More exactly, the $2-2$ state splits into the $1-1-1-1$ state when $w=w_{c}$, where

and where $x^{*}$ satisfies the implicit equation

The $1-2-1$ state, while never stable, is also tractable, with outer stack positions given by the relation

Other particular states may also have similar implicit equilibrium expressions, but their multitude makes this endeavor an impractical strategy for understanding the system for general $n$.

With five particles, we see the first case where there is multistability, by two different mechanisms. First of all, unlike $n=3$ the origin cannot stably hold a particle as we cross $w_{0}$, and the population splits into a $3-2$ state, which is necessarily asymmetric in position. Then the indifference between which stack has $3$ particles leads to bistability between two visibly different states, though they might be considered the same state up to reflection of the domain. Second, the point that the system drops to a $2-1-2$ state (on an increasing-$w$ pass) is different from the point that it jumps back to the $3-2$ state (on a decreasing-$w$ pass). We can see this in the difference between upper-left and upper-right panels in Fig. 9. This is because there is a region of bistability between $3-2$ and $2-1-2$ configurations—the loss of stability of $3-2$ happens at a higher $w$ value than the gain of stability of $2-1-2$. This hysteresis with respect to increasing and decreasing parameter $w$ is explored more generally in the “Medium-$n$ analysis” section below.

To get a sense of how the transition to the large-$n$ behavior happens, we will look at a medium-scale $n$, in particular $n=32$.

Figure 10 shows part of the bifurcation diagram for $32$ particles, which shows hints of the dynamical process by which the particles transfer between stacks. At this resolution, we can see three “connecting” branches in the top right (decreasing-$w$) figure where single particle pairs transfer from the central stack to the outer ones.

Just as with the $n=3$ and $n=4$ diagrams, there are many more equilibria than we see in Fig.  10. First of all, we only see stable equilibria due to our method of forward-time numerical integration with minutely perturbed initial conditions, so we do not see the huge number of unstable equilibria. Second, we have enforced symmetry in this simulation, so we are missing the slightly asymmetric stable states that can (and generally do) result when particles are individually free; the enforcement of symmetry is nevertheless justified as we seek a generic central pattern around which many co-stable perturbations exist. But even in the symmetric case, there is co-stability of states, which is demonstrated by the discrepancy between increasing- and decreasing-$w$ passes. Exploring these discrepancies will provide intuition about how the system behaves at higher $n$.

As Fig. 5 suggests, the overall system behavior appears to converge for large numbers of particles, under the appropriate scaling of the $w$ axis. This makes discussion of the large-$n$ limit meaningful—indeed it appears that the rapid “transitionary” bifurcations from the previous section become effectively invisible, while the “major” central stack-birth/stack-splitting bifurcations remain. There is still fractionation indifference (i.e. bands of possible stable population percentages in each stack), and the location of these major bifurcations can still vary meaningfully between forward and backward parameter continuation, as Fig. 12 shows.

We can check the stability of particular configurations like we did in Fig. 8 for $3$ particles, again using MatCont—for example, testing the stable and unstable equilibrium positions of a $129^{\textrm{th}}$ “test” particle given $128$ particles in two stacks of $64$—this is shown in Fig. 13.

In Fig. 13 we see that the test particle can align with either of the large stacks (desymmetrizing their locations imperceptibly). But we also see that the birth of stability at the origin is in fact due to a second pitchfork bifurcation with very short-lived asymmetric unstable branches which cross the outer stack positions and become stable, roughly corresponding to sitting in the trough outside the two large stacks (and in fact approaching that well location, $\pm 3s/2$, as $w\to\infty$).

For particle $j$ to be at equilibrium:

If the right hand side (RHS) is positive, particle $j$ will move right, and if negative, left. If it is zero for all particles, the system is at equilibrium. There are many particular states, such as all the particles at the origin $(0,0,\ldots,0)$, two symmetric stacks $(x^{*},x^{*},\ldots,x^{*},-x^{*},-x^{*},\ldots,-x^{*})$, etc., which may satisfy this equilibrium condition for various parameter values.

We can analytically examine the stability of the fully-stacked state at the origin, recovering the critical bifurcation value $w_{0}$ below which that state is stable—and in fact appears to be the only equilibrium.

The elements of the $j^{\textrm{th}}$ row of the Jacobian matrix $J$ for the system (9) are

which, at the origin $x_{i}=x_{j}=0$ (corresponding to the fully-stacked state), become

Due to the symmetry, we can identify all the eigenvectors:

1. 1.

   The eigenvector $v_{1}=(1,1,\ldots,1)$ corresponding to the full stack drifting left or right from the origin has eigenvalue $\lambda_{1}=-2/w^{2}$, which is always negative, indicating that the 1-stack system is stable to these types of perturbations (unsurprising based on intuition for a single particle).

2. 2.

   The other $n-1$ eigenvectors consist solely of symmetric two-particle divergence; i.e., vectors of the form $(-1,1,0,0,\ldots,0)$ with the positive $1$ in each of the other $n-1$ positions. These vectors all have eigenvalue $\lambda=-2/w^{2}+2nk^{2}/(k-1)s^{2}$. The 1-stack state is thus stable to these types of perturbations for

   $\displaystyle\frac{w}{s}$

   $\displaystyle<\frac{1}{k}\sqrt{\frac{k-1}{n}}\;.$

   (11)

   Solving for $w$, this gives us the critical parameter value $w_{0}$ for general $n$:

   $$w_{0}=\frac{s}{k}\sqrt{\frac{k-1}{n}}\;.$$

   (12)

   We note that this agrees with the $n=2$ particular case, given above after Eq. (8). This also suggests the appropriate way to rescale the $w$ axis as $n$ varies (as we have for all figures), since the structure appears to depend only on the ratio $w/w_{0}$ (at least asymptotically for large $n$).

The simplest nontrivial arbitrary-$n$ case, the two-stack state—where the population is split into two symmetric stacks of $n/2$ particles each—is quite relevant to examine since it is the dominant behavior for approximately $w_{0}<w<2w_{0}$. It is exactly solvable via a simple tweak of the logic which led us to Eq. (5), with the influence from the other “particle” being multiplied by $n/2$ while the background contribution is unchanged:

The solution can be written explicitly in terms of the Lambert W function (as done above for the special case $n=2$), now using $w_{0}$ to further simplify:

For $|x|\ll 1$, this is approximately

exactly as we saw before (all scaling is accounted for by $w_{0}$).

We seek an understanding of the second “major” bifurcation: the birth of the three-stack state near $w=2w_{0}$, with a small central stack between the two large symmetric stacks. We can easily check the curvature of the potential landscape at the origin between two equal stacks, using this as a test to identify when a particle would stably rest there. We note that this spot sees an identical contribution from all particles at $\pm x^{*}$:

so the birth of stability happens when this curvature crosses $0$, at

When we combine this condition with the two-stack equilibrium relation for $w^{2}$, Eq. (13), we get

which defies closed-form solution but which we may numerically approximate for our default parameters $k=2$ and $s=1$, yielding $x_{c}^{*}\approx 0.324$. We can see that this agrees empirically with the stack width coincident with the birth of stability in our large-$n$ figures like Figs. 4 and 12.

This approximation for $x^{*}$ in turn allows us to approximate $w_{c}$, by using either relation again. Using Eq. (13) with $k=2,s=1$, we have

This is suspiciously close to $2w_{0}$, but these approximations were done using $16$-digit precision, so it appears to indeed be distinct. We note that this value depends on $k$ (though $s$ may be scaled out as always); for example, for $k=100$ we have $w_{c,100}\approx 1.791w_{0}$.

After the central stack’s creation, the exchange of particles between central and outer stacks is complicated, since the fraction of the population at the origin influences the position of the outer stacks, which in turn influences the stable fraction at the origin. We observe empirically that the outer stacks drift apart in a roughly linear manner, which may be helpful, though we leave this exploration for future work.

We note that the increasing-$w$ sweep only sees particles at the center when they are kicked out of the outer stacks to populate it. On the decreasing-$w$ pass, however, the center hosts a larger stable population at each $w$ value reached, since continuation from the right causes the accumulation of all nearby particles at the center, where they stay until they are ejected to the outer stacks. It is perhaps counterintuitive that these particles “climb” the global potential as it narrows, but it is nevertheless true; the narrowing background potential pushes the outer stacks inward enough that the center becomes less stable, at which point it is “overstuffed” and repels some of its former constituents to join the other stacks.

The decreasing pass thus acts as a lower bound for the maximal stable fraction at the center, and the increasing pass acts as an upper bound on the minimal stable fractionation. We expect a continuous band of stable fractionations between those values, as indicated by Fig. 14. We leave further exploration of the bands of stability in the $n\to\infty$ limit for future work.

The parameter $k$ distorts the equilibrium diagram; in particular, Fig. 15 shows what happens when $k$ grows. We see that large $k$ coincides with a narrower population (smaller range on vertical axis); this is perhaps unintuitive since large $k$ makes the Ricker wavelet’s trough shallower, which intuitively means “less stability” overall, and perhaps the dominance of repulsion. However, that would imply a wider population, while in reality we see a narrower one. We might instead take away the competing explanation that the Ricker wavelets’ weaker attraction to their preferred distance means the particles don’t “buoy” their partners away from the origin as strongly, which leads to more dominance of the background potential-well containment overall.

Asymptotic analysis of Eqs.15 and 16 with large $k$ yields

which seems compatible with the large-k diagrams shown in Fig. 15.