Transport and nonequilibrium phase transitions in polygonal urn models

The study of particle transport in confined geometries, in the presence of obstacles and disorder, or of coupling with external agents concerns a large fraction of statistical mechanics literature. It has led, for instance, to the understanding of anomalous behaviors.^1–13 A classic example is afforded by single-file diffusion that concerns particles moving in narrow channels in which erratic motion is limited.^14,15 Moving defects or obstacles have also been observed to sensibly affect the residence time of a particle in 1D random walks on a lattice.¹⁶ A number of variants of the Lorentz gas, such as billiards with periodically distributed flower-shaped obstacles, have been devised to point out possible peculiarities of non-interacting particles transport, including the fractal dependence of transport coefficients on systems parameters.¹⁷ A related research line shed light on uphill currents (non-Fickian currents moving along the density gradient) in boundary driven stochastic lattice models. Uphill currents were first observed numerically, and also studied theoretically, in a 1D cellular automaton equipped with a long-range Kac potential.^18–20 The hydrodynamic limit of the model was shown to correspond to a gradient flow relative to a non-local “mesoscopic” free energy functional, which generalizes the classical Ginzburg–Landau functional. The analysis was next extended to a 2D Ising model coupled to external reservoirs, where the updating rule at the boundaries is supposed to break the condition of detailed balance. The result of such nonequilibrium driving is a magnetization (uphill) current flowing from the reservoir with lower magnetization to the reservoir with higher magnetization density.^21,22 In these models, stationary uphill currents are typically sustained by first-order phase transitions, measured by an abrupt change of the order parameter, when the temperature and the magnetization of the reservoirs are carefully chosen. Similar findings have also been discussed in the context of Zero Range processes.^23,24

Many other works have been devoted to transport problems, but an exhaustive account is beyond the scope of our paper. Here, we focus on feedback-controlled billiard dynamics that recently have been considered also under the point of view of finite size effects²⁵ and of nonequilibrium phase transitions.²⁶ While the present study focuses on polygonal urns, in order to cast our results in a proper perspective, we start by reviewing the theory developed for billiards with circular urns. In particular, a 2D billiard table consisting of two circular urns connected by a rectangular channel was taken as the container of $N$ point particles, moving on straight lines with unit speed, and undergoing elastic collisions with its walls. This amounts to an ergodic system (in fact enjoying even stronger²⁷ dynamical properties). That table was then made “active” introducing a kind of Maxwell demon in the channel that inverts the horizontal component of the particles velocities, when their number in the channel exceeds a given threshold $Θ$⁠. Such a bounce-back mechanism produces an indirect interaction between the particles that alters the standard billiard dynamics, but preserves phase space volumes. Therefore, it can be considered formally non-dissipative. The result is a sort of rectifier, which favors a nonequilibrium phase transition with non-homogeneous stationary spatial densities. The addition of a passive channel, i.e., one without bounce-back mechanisms, associates the nonequilibrium phase transition with a stationary current that takes particles from the more populated urn to the other. This discharge through the passive channel establishes a stationary flow that mimics the current in a circuit generated by a battery, because in the active channel particles flow from the less to the more populated urn.²⁸ 

Such models have been approached in probabilistic terms, relying on the ergodicity of circular urns connected by straight channels, that holds in the absence of the bounce-back mechanism.²⁷ Moreover, to make the bounce-back perturbation minimal, and the finite size effects negligible, the large $N$ and the small channels width limit were assumed. One could then rely on a uniform distribution of particles in the billiard table, as well as on a uniform distribution of the angles of their velocities. Consistently, the earlier investigation²⁶ for equal circular urns connected by a single passive channel, or, equivalently, by an active channel with $Θ≥N$⁠, showed that each particle spends, on average, the same fraction of time in each urn. Thus, a large ensemble of particles distribute themselves uniformly in such billiard tables, apart from fleeting evanescent fluctuations, that are less and less frequent as $N$ gets larger and larger. Differently, for $Θ<N$⁠, ergodicity in the absence of bounce-back implies that the number of particles in one gate will sooner or later exceed the threshold, activating the bounce-back mechanism. Such a feedback mechanism, which is more effective for smaller $Θ$ and less effective for thinner channels, promotes the transition toward steady states with non-uniform mass distribution. Indeed, increasing the number of particles in one circular urn makes easier for them to crowd its gate, while particles still come from the other urn. In fact, for $Θ$ less than a critical value (which depends on the geometrical parameters of the model), a sharp transition was observed, with a large majority of particles inside one circular urn only. Such an abrupt breaking of the left-right symmetry amounts to a nonequilibrium (first order-like) phase transition, characterized by a stationary density jump between the urns.

A difference arises between the case without passive channel and that with passive channel of finite width $wp≠0$⁠. In the first case, one circular urn gets depleted of particles, while the other urn gets more and more populated, until a steady state is reached in which the flows of particles per unit time, exchanged between the two circular urns, equalize.²⁶ In the second case, particles are able to flow between the two urns also through the passive channel, producing a net current from the more crowded to the less crowded urn.²⁸ If the ratio between channels widths $wp/wa$ is large enough, this flux balances the crowding of one urn caused by the bounce-back mechanism, keeping the homogeneous state. If $wp/wa$ is not large enough, an inhomogeneous state is realized.

To develop a theory of such phenomena, one question is whether ergodic properties may still be invoked, despite the activity of the gates. That is required for the deterministic system to be treated like a stochastic process, and the previous works proved that such an approach provides accurate descriptions of the deterministic systems. In particular, a probabilistic framework was developed for circular urns²⁶ that quantitatively yields the phase diagram of the homogeneous and inhomogeneous phases, as functions of the model parameters. The theoretical predictions, thus, obtained were confirmed by numerical simulations.

Here, we first write down the theoretical equation for the interface separating homogeneous and inhomogeneous phases in terms of a restricted set of parameters, having a clear physical interpretation. One such parameter is the threshold $Θ$⁠, responsible for the bounce-back mechanism; another parameter, called $λ$ and encoding the geometry of the system, can instead be interpreted as the typical number of particles entering a gate from the adjacent urn in the interval of time in which a particle crosses the gate. A third parameter that also appears in the interface equation corresponds to the ratio of the channel widths. We make no reference to the shape of the urns, but we trust that positions and velocities are uniformly distributed in their respective spaces. This may be imposed by hand or be a consequence of good ergodic properties implied by the billiard table.

Therefore, we here investigate whether ergodicity of the dynamics in the absence of the bounce-back mechanism is required for the phase transition to occur, considering a billiard model constituted by polygonal urns whose walls make with each other angles that are irrational with respect to $π$⁠, as well as urns whose sides make rational angles with each other. In fact, the family of rational polygons is not ergodic: given the initial one, they generate only a finite number of velocities, and we analytically evaluated their number as well as their directions. The lack of ergodicity implies that the velocity distribution is anisotropic, which makes our theory not immediately applicable. Moreover, the finiteness of the number of velocities allows the identification of periodic orbits, hence of particles trapped forever in one urn. However, periodic orbits constitute a lower dimensional subset of the phase space and are irrelevant in simulations with random initial conditions. On the contrary, irrational polygons are expected to be ergodic. A general proof is lacking,^29–31 but in our case, a standard argument proves that the velocities generated from typical initial conditions are infinitely many and dense in the unit circle.

Our numerical investigation of a billiard model with polygonal urns will ultimately confirm that, while a nonequilibrium phase transition may take place in rational polygonal billiards, it cannot be quantitatively captured by our theoretical framework, owing to the lack of ergodicity of the dynamics in the two connected urns. Thus, the picture which emerges is that: while the microscopic origin of the phase transition can be traced back to the velocity-reversal mechanism acting in the gates, a mathematical characterization of the interface region, based on a probabilistic reasoning, can be achieved only if some sort of ergodicity holds. In fact, a related question is whether the bounce-back mechanism hinders or favors ergodicity. The present study of rational polygonal billiards will reveal that the number of velocities one particle can experience is doubled by the velocity-reversal mechanism in the gates. Consequently, rational tables remain non-ergodic, even when the rate of activation of the bounce-back rule is high. We also show analytically that the velocity space is densely explored for our irrational polygons, in accord with Gutkin’s conjecture.^30,31

This ultimately clarifies the nature of the bounce-back mechanism: while it is an essential ingredient in the emergence of the phase transition, it does not alter significantly the ergodic properties of the microscopic dynamics, as originally assumed. Then, the detected major discrepancies between our theory and the behavior of rational polygons, while there is perfect accord in the case of circular urns or irrational polygons, confirms that the probabilistic reasoning requires good ergodic properties. It is, thus, interesting to note that while transport phenomena in polygonal billiards have already been studied,^32–37 the presence of the bounce-back mechanism constitutes a interesting novelty from the dynamical systems perspective.

Thanks to these facts, our investigation of phase transitions distinguishes between rational and irrational polygons. Furthermore, it accounts for the presence of the passive channel, see Fig. 1, which allows particles to flow freely, contrasting the trend toward inhomogeneous states.

This work is organized as follows: in Sec. II, we describe the model by introducing the microscopic dynamics and the observables of interest. In Sec. III, we develop the theoretical framework to characterize the phase diagram of the model. In Sec. IV, the results of the numerical simulations are reported and compared with the theoretical predictions. Conclusions are drawn in Sec. V.

Consider a billiard table $Ω$⁠, made of two identical rhomboid-shaped urns, referred to as urn 1 and urn 2 for the left and the right urns, respectively. Each urn $U⊂Ω$ is connected to two rectangular channels with widths $wa$⁠, $wp$ and lengths $ℓa$⁠, $ℓp$⁠, called active and passive channels. The openings of the channels are located in correspondence of the horizontal vertices of the two urns, while the longitudinal side of each channel extends along the $x$-axis, see Fig. 1. We also denote by $ℓr$ the length of one side of the urn and by $α$ the amplitude of the upper vertex angle. The active channel is constituted by two halves (the blue-shaded regions in Fig. 1), each of length $ℓa/2$⁠, called gates. The passive channel is also split into two halves, one located on the left of urn 1 and the other on the right of urn 2. Periodic boundary conditions are enforced along the $x$-direction at the leftmost and rightmost end points of the passive channel, thus forming a closed loop. $N$ independent point particles move in straight lines with speed $v=1$ inside $Ω$⁠. The particles collide elastically with the boundaries $∂Ω$ of the table, where they undergo specular reflections.

Depending on the amplitude of the vertex angle $α$⁠, the billiard table is classified as a rational polygon, if $α/π=p/q$ with $p,q∈N$⁠, or an irrational polygon, if $α$ is not a rational multiple of $π$⁠. Such distinction is of paramount importance in establishing the ergodic properties of the dynamics in polygonal billiards which to date is, in general, an open problem.^30,31

What distinguishes our model from a standard billiard is a dynamical rule, characterizing the dynamics in the two gates of the active channel, originally introduced in Ref. 26: the bounce-back mechanism. Such rule works as follows: when the number of particles in one gate, which are moving toward the other gate, exceeds a fixed threshold $Θ$⁠, the horizontal component of their velocity is reversed. Particles in a gate that are coming from the other gate are not affected by the mechanism and may continue their motion unaltered.

A quantitative study of the phase transition in our model starts with the definition of a suitable order parameter. To this aim, we introduce the observable called mass spread, defined as

where $N1(t)$⁠, $N2(t)$ denote the (instantaneous) number of particles residing in urn 1 and urn 2, respectively, at time $t$⁠. We then compute the absolute value of the time averaged mass spread, $χ¯$⁠, viz.,

with $T$ large enough to ensure that the time average attains an asymptotic value. The observed steady states are, thus, described by different values of $χ¯$⁠, while homogeneous states correspond to $χ¯=0$⁠, $χ¯$ close to $1$ refers to unequal distribution of particles in the urns.

Another important observable is the particle current, which is defined as follows. Let $x=0$ be the position of the vertical line separating the two gates of the active channel and consider a long time interval $[0,T]$⁠. Denote by $Δn(T)$ the number of particles that in the time $T≫1$ cross $x=0$ in the direction from urn 1 to urn 2 minus those moving in the opposite direction. Then, the net stationary current per particle is defined by

Before proceeding further with the theoretical description of the model, in the remainder of this section, we shall focus on some ergodic properties of the billiard models considered in this work. To begin with, we shall restrict our discussion to closed rhomboid-shaped polygons, i.e., urns with $wa=wp=0$⁠. Let $θ0$ be the direction of the velocity of a particle at time $0$ with respect to the positive direction of the horizontal axis and $θk$ be the corresponding direction after $k$ collisions. The sequence of velocities fully specified by $(θk)k∈N$ is obtained from the mapping $Mα:[0,2π]×{−1,+1}→[0,2π]$ defined as

where the sign $σk∈{−1,+1}$ is determined by the position of the reflection point on the surface of the urn: $σk=+1$ for the opposing top right and bottom left sides of the urn, and $σk=−1$ for the remaining two sides. It is important to apply the modulus in Eq. (4) at each iteration and not only once after $k$ iterations. A natural question, thus, regards the characterization of the allowed velocity directions in the urn. It turns out that for irrational polygons (⁠$α/π$ is an irrational number), $θ$ densely covers the whole set $[0,2π]$⁠. This is immediately seen from Eq. (4): the map $Mα$ generates an irrational rotation on the circle $R/(2πZ)$⁠; hence, orbits are dense in $[0,2π]$ and the dynamics are ergodic. This leads to a uniform exploration of the billiard table, meaning that sets of equal area are visited with equal frequency. However, the rate at which such a uniform distribution is reached sensibly depends on $θ0$⁠, cf. Fig. 2 for the rate at which the interval $[0,2π]$ is explored starting from different initial velocity angle $θ0$⁠.

Instead, for rational polygons (i.e., if $α/π=p/q$⁠, with $p,q$ some integer numbers), the angle $θ$ takes up just a finite number of values. In fact, for a given $θ0$⁠, there are at most $2q$ admissible velocity directions, and exactly $2q$ different directions if $θ0$ is irrational. Specifically, if $q$ is odd,

holds, whereas taking even $q$ in Eq. (4) yields

For instance, for a closed rational urn with $θ0=π/3$⁠, which is irrational, Eq. (5) evaluates to

for $α=π/3$ (odd $q$⁠), while Eq. (6) produces

for $α=π/2$ (even $q$⁠), which is also shown by our numerical simulations (see Figs. 6 and 13).

In the presence of the active channel (⁠$wa>0$⁠), a velocity reversal in $x$-direction corresponds to $σk=0$ in Eq. (4). This may lead to additional members in the set of possible $θ$ values. For the case of Eq. (7), the additional six angles are given by

while for the case of Eq. (8), there are no additional angles in agreement with numerical results to be presented below. Thus, as expected, dynamics in single rational urns is found to be non-ergodic. Note that for rational polygons with rational $θ0$ the above general expressions still apply, but the number of different $θ$ is eventually reduced; repeated entries have to be ignored from the list of size $2q$⁠. For example, for rational $θ0=α/2$⁠, every entry in the list appears twice, and there are only $q$ different $θ$’s.

In Sec. III, we develop a probabilistic theory for the billiard table shown in Fig. 1, which relies on the ergodic properties of the dynamics.

In this section, following ideas developed earlier,^26,28 we outline briefly the theoretical framework describing the phase transition in the system. A key ingredient in our derivation is the single-particle invariant probability density $ρ(x,v)$ to find one particle at $x∈U$ with velocity $v=(vcos⁡θ,vsin⁡θ)$⁠, where the speed $v>0$ is fixed. Analogously to Refs. 26 and 28, our derivation relies on the ergodicity of the irrational polygonal billiard constituted by two urns connected by a channel with straight walls and on the assumption that such ergodic properties are not substantially altered by the introduction of the bounce-back mechanism. This is of course favored by the smallness of the channel widths. Furthermore, given that the dynamics preserves phase space volumes, we assume that $ρ(x,v)$ is uniform in $U$ and on the interval $[0,2π]$⁠. It can, thus, be written as $ρ=(2πA)−1$⁠, with

the area of a single rhomboid-like urn. While the assumption about the uniformity of the density seems appropriate for irrational polygons, it is unjustified for rational polygons, in which only few velocity directions are in fact admitted. We may, hence, expect our theory to capture the phase transition more appropriately with irrational than with rational polygons: this will indeed be confirmed by our numerical investigation, in Sec. IV. Let us call

the stationary number of particles in the $i$th urn (⁠$T≫1$⁠), as in Eq. (2) and $τa$ (resp. $τp$⁠) is the typical time needed for a particle to cross one gate in the active (resp. passive) channel.

Owing to the uniformity of the density $ρ(x,v)$ (in the simulations, this requires large $N$⁠, otherwise unoccupied regions that defeat uniformity are necessarily too large), expressions can be derived for the stationary current flowing through the active and passive channels (details are given in the  Appendix). They are determined by $λi,a$ (or $λi,p$⁠) that can be interpreted as the typical number of particles which, in the time interval $τa$ (or $τp$⁠), enter the active (or passive) gate from the adjacent $i$th urn. We find

for $i∈{1,2}$⁠, and

For the current flowing through the passive channel, we obtain in terms of these quantities

The current flowing through the active channel is instead reduced by the action of the bounce-back mechanism in the gates and attains the modified form

where $Γ(y,x)$ is the Euler upper incomplete gamma function defined in Eq. (A6) and $Γ(y)=Γ(y,0)$⁠.

The current $J¯a$ is conventionally taken as positive if the flow goes from the left to the right urn through the active channel. Analogously, $J¯p$ is conventionally taken as positive if the flow goes from the right to the left urn through the passive channel (periodic boundary conditions are imposed).

The condition of stationarity requires the equality of the currents $J¯a$ and $J¯p$ in Eqs. (14) and (15). Therefore, letting $η:=Nτa(J¯a−J¯p)$⁠, one has

Using the definitions (12), one has $λi,p=λi,awpℓp/waℓa$⁠, and noting that $τa/τp=ℓa/ℓp$⁠, Eq. (16) can, thus, be written as

Note that here the ratio $wp/wa$ can be taken as a single parameter. That is also possible because our results can be scaled and we choose a length unit so that $wa=1$⁠. However, in the following, we continue writing $wp/wa$ to highlight its meaning as a ratio of widths.

When most of the particles reside inside the urns, e.g., in the case $N1,N2≫Θ$⁠, we have $N≈N¯1+N¯2$⁠, and the typical number of particles entering the active channel from both urns, in a time $τa/2$⁠,

obeys $λ=(λ1,a+λ2,a)/2$⁠. Consequently, $λ2,a=2λ−λ1,a$⁠, and the stationary mass spread is given by

The existence of multiple solutions of Eq. (17) signals the onset of a phase transition. One solution, as it can be easily verified by direct inspection, is $λ2,a=λ$⁠, or equivalently, $N¯2=N/2$⁠, which corresponds to the equilibrium state $χ¯=0$⁠. More solutions may appear, though, which bring the stationary dynamics out of equilibrium. Indeed, as shown in Figs. 3 and 4, for certain parameters values, $η$ changes sign in intervals not containing $λ2,a=λ$⁠. The roots in these intervals correspond to nonequilibrium states, $χ¯>0$⁠. Given the differentiability of $η$⁠, one may ask which of the steady states is stable, when $η$ vanishes more than once for $λ2,a∈[λ,2λ]$⁠.

The linear stability analysis states that a stationary state is stable if $(∂η/∂λ2,a)≤0$⁠, otherwise it is unstable.^26,28 Hence, the theoretical transition line between homogeneous and inhomogeneous phases is identified by the locus of points at which

Using Eqs. (17) and (18), the condition (20) becomes

or, equivalently,

Equation (22) is a central result of our theoretical derivation, as it yields the behavior of the system determined by the parameters $Θ$⁠, $λ$⁠, and $wp/wa$⁠. It holds not only for polygonal urns but also for circular urns, if $A$ is the area of a circular urn. Note that most (all for $wp=0$⁠) geometric parameters determining the shape of the billiard table, $ℓr$⁠, $ℓa$⁠, $α$⁠, $wa$⁠, and $wp$⁠, directly contribute to the value of $λ$⁠. Then, fixing the ratio $wp/wa$ in Eq. (22), one obtains the equation $λ=λ(Θ)$ for the line delimiting the different phases.

For $wp/wa=0$⁠, Eq. (22) has a single solution separating the only two stationary states that exist when the passive channel is absent: the homogeneous state (⁠$χ¯=0$⁠) and the totally inhomogeneous state (⁠$χ¯=1$⁠).

For positive but sufficiently small $wp/wa$⁠, the line $λ=λ(Θ)$ turns concave with a maximum at $λ=Θ+1$⁠. Then, Eq. (22) has two solutions in the interval $[λ,2λ]$ that represent the two interfaces $λ+(Θ),λ−(Θ)$⁠, which split the parameter space into three regions. The first, also occurring in the single channel model, corresponds to the homogeneous state in which the bounce-back is in-effective. The second corresponds to an inhomogeneous state, $0<χ¯<1$⁠, and a steady current in the circuit. In this case, a net current from the more to the less populated urn flows through the passive channel, while the active channel plays the role of an efm, that maintains a kind of potential difference, allowing more particles to go from the less populated urn to the more populated urn than vice versa. The third region, characterized by $χ¯=0$⁠, corresponds to a situation in which the particles exiting the urns through the passive channel are not efficiently replaced by those bouncing back in the active channel. The state is homogeneous and no net current is produced.

It is finally interesting to note that a stable homogeneous and a stable inhomogeneous state may coexist, as shown in Fig. 3(a), where two inhomogeneous stationary states appear: one is unstable and the other one is stable. The theoretically predicted behavior of $χ¯$ is reported in the plots in Fig. 5, whereas numerical evidence of the coexistence of homogeneous and inhomogeneous states will be shown in Fig. 9.

In Sec. IV, we will compare our theory with the results of the numerical simulations of the dynamics. The comparison will enable us to reveal the role of the ergodicity of the dynamics in our theory of nonequilibrium phase transitions since that depends on the geometry.

The theoretical analysis of the stationary behavior of the model, developed in Sec. III, is complemented here with a detailed computational study of the dynamics inside the table $Ω$⁠, that does not rely on the assumptions of the theory. A bounce-back event-driven algorithm has been implemented to evolve the particle trajectories basically independently, at asynchrone time steps, to allow for maximum execution speed. Each particle propagates directly, from its current position, to the time and position where it enters or exits a gate of the active channel. Not only is the period of time, thus, a variable for each particle, but also the number of its collisions with the system walls until it reaches or exists a gate (a much more time-consuming but simpler implementation runs at identical time step for all particles, at a variable time step that corresponds to the shortest time until next collision; less efficient implementations exist as well). Entrance and exit times for each single particle are recorded and analyzed, whenever the number of detected particles contemporary spending time in one of the gates exceeds the threshold $Θ$⁠, their trajectory is rewinded and the horizontal component of the velocity of those occupying the gate is reversed. This way the bounce-back mechanism is implemented exactly. A steady state detection procedure is implemented, based on the constantly recorded value of $χ(t)$⁠; in particular, the beginning of the averaging period is defined as the time instant $t$ at which the distance $||χ(t)|−|χ(0)||$ has reached its all time maximum value. Because $χ(t)$ tends to approach a stationary value within fluctuations, the whole averaging period is considered to represent the steady state. Once the steady state is detected, the simulation evolves the dynamics for an adjustable amount $T$ of time, according to the desired accuracy of the averaging procedures.

In each set of simulations, the number of simulated particles $N$ is fixed, while the geometrical features of the billiard table (i.e., $ℓr$ and $α$ for the urns, $wa$⁠, $wp$⁠, $ℓa$⁠, and $ℓp$ for the active and passive channels) and the threshold value $Θ$ are adjusted to explore the parameter space. We are going to present results for the initial state, where all particles reside in one of the urns, $χ(0)=1$⁠. This choice is computationally advantageous, as the regime of large $λ$ and small $Θ$⁠, for which $χ¯=1$ is the usual scenario, is otherwise only reached at an enormous computational effort. We have paradigmatically verified that results are insensitive to the choice of $χ(0)$⁠.

Note that the system parameters can be limited to $λ,Θ$⁠, and $wp/wa$⁠, since they determine the stationary properties of the dynamics, via Eq. (17). Thus, the assumption that $N1+N2≈N$ in the limit of small channels, and Eq. (12), imply that the length $ℓp$ of the passive channel does not belong to the set of parameters to be explored. Nevertheless, we remark that the magnitude of $ℓp$ may affect the simulation times and must also be properly tuned in order to meet the assumptions introduced by the theory. The results are presented in two different subsections: the first shows simulation results for irrational polygon with vertex angle $α=π/8$⁠, while the second refers to a rational polygonal table with $α=π/3$⁠. For both setups, two different initial distributions of particles velocities are tested, corresponding to random uniformly distributed initial velocities and to unidirectional initial velocities pointing along the direction $θ0$⁠. Sample trajectories starting from $θ0=π/3$⁠, for the two setups with irrational and rational polygonal urns, are illustrated in panels (a) and (b) of Fig. 6. Panel (d) of the same figure shows that, in the case of rational polygonal urn with $α=π/3$⁠, the values of the velocity directions after the $k$th collision on the walls, identified by the angle $θk$⁠, match the finite set of values predicted by Eq. (7). Instead, in the case of an irrational polygonal urn with $α=π/8$⁠, the values of $θk$ tend to fill up densely the interval $[0,2π]$⁠, as evidenced in panel (c).

Here, we report the results from simulations of the irrational polygon described above. Figure 7 shows the steady state mass spread $χ¯$ as a function of the parameters $λ$ and $Θ$⁠, for different values of the ratio $wp/wa$⁠.

As the plots show for both types of initial conditions, the case with a single active channel (i.e., $wp=0$⁠) behaves as theoretically expected: the inhomogeneous steady state with $χ¯=1$ is confined in the parameters space region predicted by the theory. The correspondence with the solution of the interface equation is perfect, which validates also the ergodic hypothesis underlying our reasoning. The same cannot be directly said about the configurations with $wp≠0$⁠. For instance, taking $wp/wa=0.05$⁠, none of the initial conditions of velocities seems to precisely drive the system toward the expected behavior. The agreement between theory and simulations, cf. the lower green transition line in the bottom panels of Fig. 7, is better realized at small values of $λ$⁠, namely, when the rhomboidal geometry of the urns is not significantly affected by the presence of the channels. However, also this is in accord with our theoretical derivations, which assume the channels negligible in size, so the standard billiard dynamics be only weakly altered.

The dependence of $χ¯$ on $N$⁠, when the geometrical parameters are fixed, is shown in panel (a) of Fig. 8. The plot gives $χ¯$ as a function $Θ/λ$ for different values of $N$⁠: an increase of $N$ makes the transition region sharper and sharper, thus resembling the classical behavior observed in first-order phase transitions in the large $N$ limit. In fact, this limit is necessary also for the stability of the theoretically predicted inhomogeneous stationary states, because, at finite $N$⁠, only the homogeneous state is stable. The reason is that the homogeneous state occupies by far the largest fraction of the phase space and ergodicity implies that sooner or later it is achieved. However, fluctuations away from equilibrium are possible at finite $N$ and last longer and longer as $N$ grows. This does not contradict Figs. 3 and 4, since they are based on a probabilistic theory that assumes a uniform distribution in positions and velocities. Hence, it applies only to systems made of many particles and turns more accurate for larger $N$⁠, if particles are randomly distributed. Indeed, Fig. 9 shows how inhomogeneous states considered stable by the theory turn from unstable to stable as $N$ grows, and within the range of $N$ values in which they exist (changing $N$ with other parameters fixed may change the number of steady states), cf. Fig. 4. The same holds for circular urns.

Letting $A=2A+ℓawa+ℓpwp$ be the total area of the billiard table and $N/A$ be the mean particle number density, Figs. 10 and 11 display the reduced stationary particle density profiles, i.e., local number density multiplied by $A/N$⁠, and the velocity profiles in the urns and in the channels, in presence [panels (a), (c), (e)], and in absence [panels (b), (d), (f)] of the bounce-back mechanism. For different values of the parameters $(Θ,λ)$ at $wp/wa=0$ and at $wp/wa=0.05$⁠, a comparison of panels (a) and (b), as well as of panels (c) and (d) of Fig. 10, shows that a uniform distribution of particles is always obtained in the absence of bounce-back, but even in the presence of bounce-back events, the distribution can be uniform, if $λ$ and $Θ$ are properly tuned. Furthermore, the histograms of the velocities in the urns are uniform in $[0,2π]$⁠, as required for the application of the probabilistic method of Sec. III. As discussed in Sec. II, this owes to the fact that the map $Mα$ in Eq. (4) is ergodic for irrational polygonal urns. When the bounce-back mechanism is efficeintly at work and produces a phase transition, the distribution of velocities in the channels is naturally altered.

An analogous numerical investigation has been carried out also for irrational polygonal billiards equipped with a passive channel, with $wp/wa=0.05$⁠. In this case, a stationary current may be established in the circuit constituted by urns and channels, with the active channel acting as an emf. Our simulations show that both the density and the velocity histograms referring to the crowded urn (for which a meaningful statistics is available) remain approximately uniform, even in the presence of a stationary current, cf. Fig. 11.

Here, we report our results, obtained from simulations of rational polygons dynamics. Figure 12 shows the steady state mass spread $χ¯$ as a function of the normalized parameters $λ/N$ and $Θ/N$⁠, for different values of the ratio $wp/wa$ with $α=π/3$⁠. Left and right panels in the figure refer to unidirectional initial velocities and random initial velocities, respectively; top and bottom panels refer to tables with $wp/wa=0$ and $wp/wa=0.05$⁠, respectively. Unlike the irrational case, our rational case with $wp=0$ matches the theoretical predictions only when a uniform initial velocity distribution is simulated, cf. panel (b) in Fig. 12. Indeed, as mentioned in Sec. II the map $Mα$ in Eq. (4) is non-ergodic for rational rhomboidal urns. However, as revealed by panels (b) and (d), the difficulty for the application of our theory, which assumes that all velocities are available, is overcome if the initial distribution of velocities is random and $N$ is large. Deviations from the theoretical predictions are instead visible in panel (a) of Fig. 12 and become even more severe in panel (c), where the phase transition foreseen by the theory is not even detected. This result clarifies further the role of ergodicity and of the bounce-back mechanism. Despite being a key ingredient, the bounce-back does not suffice, per se, to produce the phase transition, for some degree of uniformity of the velocity distribution is also required.

As also observed with irrational urns, the addition of the passive channel, which may give rise to a stationary current flowing through the billiard table, leads in panel (d) of Fig. 12 to some tiny deviations from the theoretical lower interface lines at large values of $λ$ and $Θ$⁠.

Figures 13 and 14 confirm the presence of a finite number of velocities in the steady state dynamics: the position of the peaks is well described by the formula given in Eq. (7), as also seen in panel (c) of Fig. 6. An inspection of the velocity histograms in the left and the right columns of both figures also highlights the role of the bounce-back mechanism in doubling the number of velocity directions. Hence, the mechanism is not suitable to produce the uniformity of the velocity distribution. The behavior of $χ¯$ as a function of $Θ/λ$⁠, for a rational rhomboidal urn with $α=π/3$⁠, is shown in panel (b) of Fig. 8, which evidences the sharpening of the transition region when $N$ grows, with random initial velocities. Figure 15 shows the results of an analogous investigation for the peculiar rhomboidal urn with $α=π/2$ (i.e., $q$ is even). As in Fig. 12, simulations confirm that an (in fact remarkable) agreement with the theory is achieved only in the presence of random initial velocities. Furthermore, as predicted by Eq. (8), the number of velocity directions, visible in Fig. 16, is less than with $α=π/3$ (i.e., $q$ odd). Another aspect which is worth highlighting is that, with $α=π/2$⁠, the bounce-back mechanism (Fig. 17) does not lead to a doubling of the velocity directions as previously noticed in Figs. 13 and 14 with $α=π/3$⁠.

In this work, we studied the deterministic dynamics of $N$ particles moving in a polygonal billiard table. We discussed the interplay of two important ingredients of the model: the bounce-back mechanism in the gates of the active channel and the geometry of the rhomboidal urns, which sensibly affects the ergodic properties of the dynamics. Our investigation of the steady state dynamics revealed that non-equilibrium phase transitions predicted by our probabilistic framework may occur in the large $N$ limit, as in the case of circular urns.^26,28 The phase transitions are observed in the behavior of the stationary value of the observable $χ¯$⁠, called mass spread, which may sharply switch from $0$ (corresponding to a homogeneous state, in which the urns are equally crowded) to a value close to $1$ (corresponding to an inhomogeneous regime with differently populated urns). The phase transition is triggered by the bounce-back mechanism, which kicks back the particles coming from an urn whenever their number in a gate exceeds a threshold $Θ$⁠. In the absence of the aforementioned mechanism, particles are independent, and the only steady state is the one that occupies the largest fraction of phase space, i.e., the homogeneous state. To a certain extent, the perturbation introduced by the bounce-back mechanism in a gate reminds the action of a classical Maxwell demon. As this demon acts on particles coming from both urns, one cannot foresee which of the two urns will eventually be crowded, if the initial state is $χ(0)=0$ and the parameters imply $χ¯=1$⁠.

The ergodic properties of the unperturbed dynamics constitute another key ingredient in the onset of phase transitions, and in the application of the stochastic theory. Indeed, this theory predicts such transitions under the assumption that a uniform distribution of positions and velocities be realized. To make sense, this in turn requires large $N$⁠, because there can be no uniform occupation of the positions and velocities spaces with a few particles. We have addressed, first, the study of irrational rhomboid-shaped urns and proved that the corresponding dynamics is ergodic. Our simulations with large $N$ not only confirmed the presence of the phase transition, but also that it is correctly described by our theory. With rational polygons, the picture changes. Depending on the geometry of the system (which may comprise, or exclude, the presence of the passive channel), our numerical simulations detected, in some cases, a phase transition, but the agreement with the theory turned out to be non-optimal. Moreover, this required the initial velocity distribution to be uniform, to beat the lack of ergodicity in rational urns, proved in Sec. II.

This clarifies the role of the various features of the microscopic dynamics in the description of collective phenomena, such as non-equilibrium phase transitions, and paves the way to further investigations, including the study of 3D billiard dynamics in the presence of the bounce-back mechanism. It is worth stressing that, for insufficiently large $N$⁠, our model correctly displays some peculiar features dictated by the microscopic dynamics and by the initial conditions, such as the migration in time from one state to another, among those predicted to be stable by the theory.

It is also interesting to highlight the different role played by the rectangular channels in billiard tables constituted by either circular or rhomboidal urns, as regards the onset of ergodicity and stronger dynamic properties. In fact, the circular urns connected by a rectangular channel (without the bounce-back mechanism) yield chaotic dynamics (K-system), while the polygonal urns have vanishing Lyapunov exponents. In the first case, the presence of mixing corroborates the use of a stochastic description of the deterministic model, as highlighted in Ref. 26.

A relevant variant of the model discussed in this paper could be obtained introducing circular scatterers inside the rhomboid-shaped urns, thus allowing chaotic dynamics.³⁸ The onset of a phase transition in such a model, induced by the bounce-back mechanism, is an interesting problem which is left to future investigation.

Unraveling connections with biological systems also offers promising scenarios. Indeed, the bounce-back mechanism exhibits an instance of a positive feedback control which is often observed in biology. One example is provided by the phenomenon of blood clotting: an injured tissue releases signal chemicals which activate platelets in the blood; the platelets, in turn, release a larger amount of chemicals, which activate more platelets, and so on. The foregoing process, thus, results in a cascade eventually leading to the formation of a blood clot.^39–43 Another example of positive feedback loop comes from the Hodgkin cycle in membrane biology: the leakage of sodium ions through sodium channels, caused by the membrane of a nerve fiber, alters the membrane potential.^44–48 This, in turn, causes a further opening of the channels, which leads to an increased flow of sodium ions. The process gives rise to a rapid increase in sodium leakage through the channels, which finally results in the generation of the nerve signal. The polygonal billiards considered in this work constitute models amenable to detailed analysis that represent a first approximation to such important phenomena.

L.R. acknowledges partial funding by the Ministero dell’Istruzione, dell’Università e della Ricerca (MIUR), Italy under Grant No. E11G18000350001 “Dipartimenti di Eccellenza 2018–2022.” L.R. further acknowledges support from the Italian National Group of Mathematical Physics (GNFM) of INdAM.

The authors have no conflicts to disclose.

Emilio N. M. Cirillo: Conceptualization (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Software (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). Matteo Colangeli: Conceptualization (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Project administration (equal); Supervision (equal); Validation (equal); Writing – original draft (equal); Writing – review & editing (equal). Antonio Di Francesco: Conceptualization (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Software (equal); Supervision (equal); Writing – review & editing (equal). Martin Kröger: Conceptualization (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Software (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). Lamberto Rondoni: Conceptualization (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Supervision (equal); Writing – review & editing (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.

Here, we derive explicit formulas for the stationary currents in terms of geometrical and probabilistic arguments.^26,28 Let $A$ be the area of a rhomboidal urn, expressed by Eq. (10), in the case of two channels. Let $ρ(x,v)=1/2πA$ denotes the probability density for a particle with velocity $v=(vcos⁡θ,vsin⁡θ)$ to be located at the point $x$ inside an urn, where $θ∈[0,2π]$⁠, and the dimensionless speed $v>0$ is fixed. Assume that the probability density is identical in the two urns.

Consider now a generic channel of width $w$ and length $ℓ$ connected to the urn. The assumed uniformity of the single-particle probability density yields, via a geometrical argument, the typical number of particles that, in a small time interval of length $δ>0$⁠, leave the $i$th urn to enter the adjacent gate of width $w$ and length $ℓ/2$⁠. Such a number reads

where $N¯i$ is the average number of particles residing in the $i$th urn. In Eq. (A1), we identify the quantity $wvδ/πA$ with the probability $pδ$ that one particle in the urn enters the adjacent gate in the time interval $δ$⁠.

Now, consider the typical time $τ$ needed for a particle to cross a gate of length $ℓ/2$⁠, i.e.,

which is computed from the knowledge of the length of the gate and the average horizontal component of particles velocity, $⟨vx⟩=2v/π$⁠.⁴⁹ The interval $τ$ is then partitioned into smaller intervals of length $δ$⁠. The probability that a number $s≪N¯i$ of particles enter a gate from urn $i$ in the time interval of length $τ$ is

Since $pδ/δ$ does not depend on $δ$⁠, the Poisson limit theorem, for small $δ$⁠, implies that the probability (A3) is well approximated by $λise−λi/s!$ with

The variable $λi$ represents the typical number of particles that enter the gate adjacent to the $i$th urn in the time interval $τ$⁠.

The probability that at most $Θ$ particles, with $Θ$ small with respect to $N$⁠, enter this gate adjacent to the $i$th urn during the time interval of width $τ$ is

where we recall the definition of Euler’s upper incomplete gamma function

Let us now consider a larger time interval of width $t$ and let us partition it in small intervals of width $τ$⁠. Typically, only in $Pτ(i)t/τ$ intervals, out of the total $t/τ$⁠, the particles entering the gate will not be affected by the type of channel, active vs passive. Therefore, the typical number of particles that cross the gate per unit time can be estimated as the typical number of particles that enter the gate in an interval of length $τ$⁠, multiplied times $Pτ(i)t/τ$⁠, i.e.,

If, on the other hand, $Θ>N$⁠, the probability that at most $Θ$ particles enter a gate adjacent to the $i$th urn during the time interval of width $τ$ is just unity, $Pτ(i)(Θ)=1$⁠, and the typical number of particles that cross the gate in the time interval $τ$ equals the typical number that has entered the gate in the same time interval from the $i$th urn, which is $λi$⁠. Hence, the current leaving the $i$th urn into the channel reads $λi/τ$ in that case.

These general considerations immediately lead to the expressions for the specific currents for the system with active and passive channels, with widths $wa$ and $wp$⁠. The dynamics in the active channel, divided into two gates, is driven by the bounce-back mechanism (if $Θ<N$⁠), while in the passive channel the particles always get through.

Upon replacing $w$⁠, $τ$⁠, $λi$⁠, and $ℓ$ in the above expressions by their counterparts indexed by subscripts $p$ and $a$ for the passive and active channels, cf., Eqs. (12) and (13), the specific average net currents $J¯p$ and $J¯a$ in the passive and active channels, accounting for particles leaving an urn and also for those coming from the other urn, are therefore given by

and

which concludes the derivation of Eqs. (14) and (15).