Computer-assisted proof of ergodicity breaking in expanding coupled maps

TL;DR

This paper presents a computer-assisted method to rigorously demonstrate ergodicity breaking in a family of coupled expanding maps, providing evidence of phase transitions in deterministic chaotic systems.

Contribution

It introduces a novel numerical algorithm to construct asymmetric ergodic components, advancing the understanding of phase transitions in deterministic dynamical systems.

Findings

Empirical evidence of symmetry breaking with increasing coupling strength.

Algorithm successfully constructs ergodic components for small N.

Indicates phase transitions are provable for systems with arbitrary particles.

Abstract

From a dynamical viewpoint, basic phase transitions of statistical mechanics can be regarded as a breaking of ergodicity. While many random models exhibiting such transitions at the thermodynamics limit exist, finite-dimensional examples with deterministic dynamics on a chaotic attractor are rare, if at all existent. Here, the dynamics of a family of $N$ coupled expanding circle maps is investigated in a parameter regime where absolutely continuous invariant measures are known to exist. At first, empirical evidence is given of symmetry breaking of the ergodic components upon increase of the coupling strength, suggesting that breaking of ergodicity should occur for every integer $N>2$. Then, a numerical algorithm is proposed which aims to rigorously construct asymmetric ergodic components of positive Lebesgue measure. Due to the explosive growth of the required computational resources, the algorithm successfully terminates for small values of $N$ only. However, this approach shows that phase transitions should be provable for systems of arbitrary number of particles with erratic dynamics, in a purely deterministic setting, without any reference to random processes.

Figures

Tables

Table 1. Table 1: Data summary of InAsUP’s exact numerical construction

Main results
D	$N_D$	$ϵ$	$|$Cyl$|$	Succ. ratio	$\#$ InAsUP	CPU time
2	3	0.44	5	24/24	22 (43)	1.5ms (3.5ms)
3	13	0.4	9	11584/11706	1250 (9100)	0.6s (13.9s)
4	75	0.47	11	355/373∗	3.6 (6.2)$\times {10}^5$	1.4h (5.2h)
5	541	0.44	7	3/3∗	6 (10)$\times {10}^6$	120h (200h)

Table 2. Table 2: Cardinality of the sets Λ i subscript Λ 𝑖 \Lambda_{i} involved in the optimisation procedure for the collection α 𝛼 \alpha involved in ( 4.1 ).

D	2	3	4	5	6
$\#{\Lambda }_i$	2	5	16	65	326

Equations

$$(F_{N,ϵ}(u){)}_i=2u_i+\frac{2ϵ}{N}\sum _{j=1}^{N}g(u_j-u_i)\text{mod}1,1\le i\le N,$$

$$h(u)=\left\{\begin{array}[]{ccl}\lfloor u+\frac{1}{2}\rfloor&\text{if}&u\not\in\frac{1}{2}+\mathbb{Z}\\ 0&\text{if}&u\in\frac{1}{2}+\mathbb{Z}.\end{array}\right.$$

$$(D{F}_{N,ϵ}v{)}_i=2(1-ϵ)v_i+\frac{2ϵ}{N}\sum _{j=1}^N{v}_j,1\le i\le N,$$

$$(\pi_{N}u)_{i}=\left\{\begin{array}[]{ccl}u_{i}-u_{i+1}\ \text{mod}\ 1&\text{if}&1\leq i\leq N-1\\ \sum_{j=1}^{N}u_{j}\ \text{mod}\ 1&\text{if}&i=N\end{array}\right.$$

$$(B_D(x){)}_i=2h(x_i)+\sum _{j=1}^{i-1}h(\sum _{k=j}^i{x}_k)-h(\sum _{k=j}^{i-1}{x}_k)+\sum _{j=i+1}^{D}h(\sum _{k=i}^j{x}_k)-h(\sum _{k=i+1}^j{x}_k).$$

$$G_{D,ϵ}(P)\subset P\text{and}P\cap -{\text{Id}}_{{\mathbb{T}}^D}(P)=\varnothing .$$

$$\bigcup _{a=1}^{N_D}{A}_a\text{mod}0,$$

$$G_{D,ϵ}(P^t)=\bigcup _i{G}_{D,ϵ}{|}_{A_{a_i}}(P^t\cap A_{a_i})\text{mod}0$$

$$A_{a_0\dots a_{\ell }}=\bigcap _{k=0}^{\ell }{G}_{D,ϵ}^{-k}(A_{a_k}).$$

$$m={\left(({\underline{m}}_{i,j},{\overline{m}}_{i,j}{)}_{i=1}^j\right)}_{j=1}^D\in {\mathbb{R}}^{D(D+1)}\text{for which}{\underline{m}}_{i,j}<{\overline{m}}_{i,j}\forall i\le j,$$

$$P_m=\left\{x\in {\mathbb{R}}^D:{\underline{m}}_{i,j}<\sum _{k=i}^j{x}_k<{\overline{m}}_{i,j},1\le i\le j\le D\right\},$$

$$P_m\cup P_{m^{\prime }}=P_{m^{\prime \prime }}\text{mod}0.$$

$$P_m=P_{m_{\text{in}}}\cup P_{m_{\text{out}}}\text{mod}0$$

$$P_m=P_{m_{\text{in}}}\cup P_{m_{\text{out},1}}\cup P_{m_{\text{out},2}}\text{mod}0$$

$$P_m^{\alpha }=\left\{x\in {\mathbb{R}}^D:{\underline{m}}_i<\sum _{j=1}^D{\alpha }_{ij}{x}_j<{\overline{m}}_i,1\le i\le E\right\}.$$

$$\sum _{k=1}^E{\lambda }_k{\alpha }_{kj}={\alpha }_{ij},1\le j\le D.$$

$$\underline{e}_{k}(\lambda,m)=\left\{\begin{array}[]{ccl}\underline{m}_{k}&\text{if}&\lambda\geq 0\\ \overline{m}_{k}&\text{if}&\lambda<0\end{array}\right.\quad\text{and}\quad\overline{e}_{k}(\lambda,m)=\left\{\begin{array}[]{ccl}\overline{m}_{k}&\text{if}&\lambda\geq 0\\ \underline{m}_{k}&\text{if}&\lambda<0\end{array}\right.$$

$${\underline{O(m)}}_i=\underset{({\lambda }_k)\in {\Lambda }_i}{\max }\sum _{k=1}^E{\lambda }_k{\underline{e}}_k({\lambda }_k,m)\text{and}{\overline{O(m)}}_i=\underset{({\lambda }_k)\in {\Lambda }_i}{\min }\sum _{k=1}^E{\lambda }_k{\overline{e}}_k({\lambda }_k,m).$$

$$\underset{\genfrac{}{}{0ex}{}{({\lambda }_k)\in {\Lambda }_i}{({\lambda }_k)\ne ({\delta }_{k,i})}}{\max }\sum _{k=1}^E{\lambda }_k{\underline{e}}_k({\lambda }_k,O(m))<{\underline{O(m)}}_i\left(\text{resp.}{\overline{O(m)}}_i<\underset{\genfrac{}{}{0ex}{}{({\lambda }_k)\in {\Lambda }_i}{({\lambda }_k)\ne ({\delta }_{k,i})}}{\min }\sum _{k=1}^E{\lambda }_k{\overline{e}}_k({\lambda }_k,m).\right)$$

$${\overline{\underline{{\Gamma }_{D,ϵ,a}(m)}}}_{i,j}=2(1-ϵ){\overline{\underline{m}}}_{i,j}+\frac{2ϵ}{D+1}\sum _{k=i}^j{B}_{a,k},\forall 1\le i\le j\le D.$$

$$P_m=\bigcup _a{P}_{m\cap m_a}\text{mod}0,$$

$${\Gamma }_{D,ϵ}(m)=\bigcup _a{\Gamma }_{D,ϵ,a}(m\cap m_a)$$

$$O\circ {\Gamma }_{D,ϵ,a}={\Gamma }_{D,ϵ,a}\circ O.$$

$${\Gamma }_{D,ϵ,a}({\mathbb{Q}}^{D(D+1)})\subset {\mathbb{Q}}^{D(D+1)}\text{and}{\Gamma }_{D,ϵ,a}^{-1}({\mathbb{Q}}^{D(D+1)})\subset {\mathbb{Q}}^{D(D+1)},$$

$$\left\{\begin{array}[]{l}(B_{2,\epsilon}(x))_{1}=2h(x_{1})-h(x_{2})+h(x_{1}+x_{2})\\ (B_{2,\epsilon}(x))_{2}=2h(x_{2})-h(x_{1})+h(x_{1}+x_{2})\end{array}\right.$$

$$-\text{Id}x=(-x_1,-x_2)(\text{mod}1),$$

$${\sigma }_{321}x=(-x_2,-x_1)(\text{mod}1),$$

$${\sigma }_{213}x=(-x_1,x_1+x_2)(\text{mod}1)\left(\text{resp.}{\sigma }_{132}x=(x_1+x_2,-x_2)(\text{mod}1)\right),$$

$$\{$$

$$-{\sigma }_{213}\circ {\sigma }_{321},-{\sigma }_{321}\circ {\sigma }_{213}\}.$$

$$\left\{\begin{array}[]{l}(B_{3,\epsilon}(x))_{1}=2h(x_{1})-h(x_{2})+h(x_{1}+x_{2})-h(x_{2}+x_{3})+h(x_{1}+x_{2}+x_{3})\\ (B_{3,\epsilon}(x))_{2}=2h(x_{2})-h(x_{1})-h(x_{3})+h(x_{1}+x_{2})+h(x_{2}+x_{3})\\ (B_{3,\epsilon}(x))_{2}=2h(x_{3})-h(x_{2})-h(x_{1}+x_{2})+h(x_{2}+x_{3})+h(x_{1}+x_{2}+x_{3})\end{array}\right.$$

Full text

1 Introduction

In statistical mechanics, the term phase transition originally refers to the emergence of multiple Gibbs states for the Hamiltonian associated with a collective system [29], most notably accompanied with symmetry breaking, as in the Ising model. While this notion has been formulated for equilibrium measures, from the dynamical viewpoint, it conveys some breakdown of ergodicity in a related dynamical process [33]. Using out-of-equilibrium procedures, such ergodicity failures have been rigorously proved for suitably designed Markov processes and probabilistic cellular automata (PCA) [16, 23]. More precisely, Metropolis rules, Glauber dynamics and the like have been designed to account for relaxation to pre-constructed equilibrium states, including the case when several such states coexist [27].

Notwithstanding the success of this dynamical approach to phase transitions, the fact that classical mechanics is ruled by deterministic laws of motion calls for evidences in purely deterministic systems, irrespective of any consideration of random processes. Barring the use of out-of-equilibrium procedures, can chaotic attractors in deterministic analogues of random models of interacting particles systems exhibit ergodicity breaking (associated with symmetry breaking)? Despite having generated considerable attention, this question still eludes a satisfactory response that would exclude numerical flaws or unverifiable theoretical assumptions, even in basic examples such as networks of coupled expanding or hyperbolic maps.

Indeed, some examples of infinite lattices of interacting chaotic maps have been designed so that their dual dynamics acting on measures consist of PCA exhibiting phase transitions [14, 24]. Thus, out-of-equilibrium approaches can in principle be lifted to deterministic dynamical systems. However, this operation requires explicit knowledge of a coupling-intensity-independent Markov partition, in particular one that ensures a pre-selected Markov chain/PCA for the dual measure evolution. Yet, such a complete understanding of Markovian properties is rare for realistic deterministic systems. This is especially the case of chaotic collective systems with homogenizing interactions [3], which typically fail to fit the standard assumptions of the theory of dynamical systems, such as being diffeomorphisms. In general, very little is known about the symbolic dynamics, and, if at all, only for weak interaction intensity [20, 9]. This shortcoming calls for verifications of non-ergodicity in a purely deterministic setting, that are independent of any knowledge of the associated symbolic dynamics.

In this setting, various studies have reported changes in the global dynamics of coupled chaotic maps as their interaction strength is increased. Infinite lattice examples have been provided, for which the invariant densities associated with the Perron-Frobenius operator undergo bifurcations. In particular, an analogous transition to the one in the Curie-Weiss model of statistical mechanics, which only affects the dynamics at the thermodynamics limit, has been proved for the model in [6]. Moreover, convergence to a point mass for strong coupling has been established in [5] for a mean-field model acting on measures on the circle.

This convergence is reminiscent of phenomenological changes in finite systems. However, such alterations have always been observed to be preceded by a reduction in the Lyapunov dimension. Further, these changes have repeatedly resulted in stationary or periodic behaviors of the spatial averages of symmetry-related observables, see e.g. [4, 8, 21, 25]. This phenomenology is comparable with synchronization-like scenarios in which trajectories asymptotically shrink to lower dimensional subspaces. Hardly compatible with a coupling-independent symbolic grammar, they cast doubt on the nature of phase transitions in dynamical systems: do such transitions only take place at the thermodynamics limit? If not, do they necessarily require a lowering of the dimension of the attractor, or can they take place while hyperbolicity properties remain unchanged, and in absent knowledge of a Markov partition?

To address these issues, we consider a family of piecewise affine mappings $F_{N,\epsilon}$ of the $N$-dimensional torus $\mathbb{T}^{N}$, which mimic interacting particle dynamics driven by chaotic individual stirring and homogenizing interactions. Interactions are of mean-field type (all-to-all coupling) with adjustable strength $\epsilon$ (Section 2).111Details of the considerations in this paragraph and the following ones will be given below. Moreover, the mappings $F_{N,\epsilon}$ have many symmetries, most notably, they commute with flipping the sign of all coordinates. For $\epsilon=0$, the units are decoupled and evolve independently. For $\epsilon\in[0,\tfrac{1}{2})$, all mappings are expanding and must support (Lebesgue) absolutely continuous invariant measures (ACIM). Moreover, each ACIM must be ergodic - and hence invariant under the action of any symmetry - when $\epsilon$ is small enough. For such mappings, it is therefore natural to (mathematically) examine ergodicity persistence or its failure in this expanding regime, together with accompanying symmetry features.

Preliminary evidences, based on simulations of trajectories, are presented in Section 3. Completing previous findings in [12], they reveal that for each $N>2$, ergodicity is broken via symmetry breaking, for $\epsilon$ sufficiently large. For $N=3$ and 4, this phenomenology has been confirmed by analytic demonstrations [12, 30, 31]. For convenience, the proofs actually deal with $D=(N-1)$-dimensional mappings (denoted by $G_{D,\epsilon}$ below) whose ergodicity and its failure coincide with $F_{N,\epsilon}$. While these proofs seem to have captured essential information about the dynamics, their approach, which relies on inspiration from direct observations of trajectories in phase space, is hardly applicable when $D$ is large.

In order to address large $D$ (or $N$), we propose to use instead computer-assisted proof. To that goal, an algorithm is introduced (Section 4) that aims to rigorously construct asymmetric $G_{D,\epsilon}$-invariant sets containing ACIM. The algorithm relies on empirical input from simulations and is based on two important properties.222Likewise, see the extended section 4.2 and additional information in the appendices for details and numerical implementation. One characteristics is that the dynamics of suitable-for-the-proof polytopes can be encoded into one on vectors in $\mathbb{R}^{D(D+1)}$. In other words, the dynamics of suitable sets can be captured by a limited number of real variables. The other one is that computations can be designed so that, when $\epsilon$ is itself a rational number, they involve rational numbers only. Results of the algorithmic construction are given for $D$ up to 5 (Section 4.1). However, they clearly indicate that ergodicity breaking should be provable for arbitrary $D$. A discussion about possible improvements and alternative proofs is provided in Section 5.

2 Expanding systems of piecewise affine globally coupled maps

The mappings $F_{N,\epsilon}$ are defined as follows [12]

[TABLE]

where $g$ represents pairwise elastic interactions on the circle [19] and is defined by $g(u)=u-h(u)$ for all $u\in\mathbb{T}^{1}$ with

[TABLE]

Here, $\lfloor\cdot\rfloor$ is the floor function. Hence $g$ is piecewise affine, with slope 1 and discontinuities at all points of $\frac{1}{2}+\mathbb{Z}$.

Mean field coupling implies $F_{N,\epsilon}\circ\sigma=\sigma\circ F_{N,\epsilon}$ for every $\sigma\in\Pi_{N}$, where $\Pi_{N}$ is the group of permutations of $\{1,\dots,N\}$. The symmetry $g(-u)=-g(u)\ \text{mod}\ 1$ implies commutativity with the sign flip $-\text{\rm Id}_{\mathbb{T}^{N}}$.

The mappings $F_{N,\epsilon}$ are non-singular. Therefore, considerations about ergodicity of their ACIM can ignore the dynamics of sets of vanishing Lesbegue measure, and in particular, the orbits of discontinuity sets. Away from these discontinuities, $F_{N,\epsilon}$ is piecewise affine with constant derivative

[TABLE]

from which it follows that $F_{N,\epsilon}$ is expanding for $\epsilon<\frac{1}{2}$. As a consequence, its Milnor attractor [2, 26] in this regime must consist of a finite union of Lebesgue ergodic components [32], viz. the attractor of almost every trajectory must be a set of positive Lebesgue measure (thereby excluding any dimension reduction).

Focusing on this expanding regime, as mentioned above, we aim to address ergodicity of the attractor, that is whether the Lebesgue ergodic components are unique or not. Up to semi-conjugacy, this question can be examined in a more convenient family of piecewise affine mappings of the $D=(N-1)$-dimensional torus. The reduced mappings $G_{D,\epsilon}$ (defined below) have equivalent features to the originals. Namely, their symmetry group is isomorphic to $\mathbb{Z}_{2}\times\Pi_{D+1}$ and in particular, we have $G_{D,\epsilon}\circ-\text{\rm Id}_{\mathbb{T}^{D}}=-\text{\rm Id}_{\mathbb{T}^{D}}\circ G_{D,\epsilon}$. Furthermore, their asymptotic dynamics for $\epsilon\in[0,\frac{1}{2})$ must also lie in finitely many ergodic components of positive Lebesgue measure.

More precisely, the transformation $\pi_{N}$ of $\mathbb{T}^{N}$ defined by [31]

[TABLE]

(semi-)conjugates $F_{N,\epsilon}$ to the direct product $G_{N-1,\epsilon}\times F_{1,\epsilon}$ (viz. $\pi_{N}\circ F_{N,\epsilon}=(G_{N-1,\epsilon}\times F_{1,\epsilon})\circ\pi_{N}$). The mapping $F_{1,\epsilon}(u)=2u\ \text{mod}\ 1$ (which acts on the sum coordinate $(\pi_{N}u)_{N}$) does not depend on $\epsilon$ and is ergodic with respect to the Lebesgue measure on $\mathbb{T}^{1}$. Therefore, any failure of ergodicity for $F_{N,\epsilon}$ has to be concomitant with the same phenomenon for $G_{N-1,\epsilon}$.

The mapping $G_{N-1,\epsilon}$ does not involve the sum coordinate. Moreover, its (constant) derivative conveniently turns out to be a multiple of the identity in $\mathbb{T}^{N-1}$, i.e. $DG_{N-1,\epsilon}=2(1-\epsilon)\text{\rm Id}_{\mathbb{T}^{N-1}}$. Explicitly, we have $G_{D,\epsilon}=2(1-\epsilon)\text{\rm Id}_{\mathbb{T}^{D}}+\tfrac{2\epsilon}{D+1}B_{D}\ \text{mod}\ 1$, where for $i\in\{1,\dots,D\}$

[TABLE]

Perturbative arguments at the uncoupled limit $\epsilon=0$, applied to the related transfer operator acting on measure densities [17], demonstrate ergodicity for $\epsilon>0$ sufficiently small, for each integer $D$. Instead, the dynamics is not expanding at the limit $\epsilon=\tfrac{1}{2}$. All mappings $G_{D,\frac{1}{2}}$ consist of piecewise isometries and their global dynamics can be hardly described. In any case, continuation arguments appear inapplicable at this limit. In order to evaluate ergodicity or its failure at the upper end of the expanding regime, we therefore opted to collect numerical evidences.

3 Empirical results from numerical trajectories

This section presents evidences of ergodicity breaking that were obtained from numerical simulations of trajectories. The hints are of two types: rendering of trajectories when $D\leq 3$ and order parameter estimates for arbitrary $D\geq 2$.333For $D=1$, no ergodicity failure occurs since the Milnor attractor of $G_{\epsilon,1}$ is transitive, and hence ergodic, for every $\epsilon\in[0,\frac{1}{2})$ [12, 31]. We begin by presenting evidences of the first type.

In low dimension, direct visualization of asymptotic-orbit traces in phase space offers a straightforward evaluation of ergodicity/symmetry and their failure. On Fig. 1, late iterates of orbits started from representative initial conditions, and their images under symmetries, are plotted for $\epsilon$ across $[0,\frac{1}{2})$. Initial conditions were selected using random sampling of phase space, in a way to render, if not all, then the most essential attractor features, with an emphasis on detecting asymmetry.

Although the bifurcation scenarios differ for the two cases, the pictures reveal that when increasing $\epsilon$ from 0, ergodicity persists for $\epsilon$ up to some $\epsilon_{D}$ and then fails beyond that threshold. In the $D=2$ case, the transitive and fully symmetric attractor continuously splits at $\epsilon=\epsilon_{2}$ into 6 disjoint and asymmetric invariant pieces. Each emerging piece breaks all map symmetries except one (Appendix A). In the $D=3$ case, a fully symmetric invariant set exists throughout the expanding domain. In addition, 6 asymmetric invariant components discontinuously appear at $\epsilon=\epsilon_{3}$, away from the symmetric set, and persist from thereon as $\epsilon$ continues to increase. Then, at $\epsilon\simeq 0.437$, this group of partly asymmetric orbits is augmented in a similarly discontinuous way by an additional analogously persisting group, composed of 8 asymmetric orbits (see the involved symmetries in Appendix B). Despite that the phenomenology differs in the two cases, the asymmetric components always appear to be disjoint from their image under the sign-flip, ie. the $\mathbb{Z}_{2}$-symmetry generated by $-\text{Id}_{\mathbb{T}^{D}}$ is systematically broken.

The analytic proofs of ergodicity breaking in [12, 30, 31] established the existence of so-called InAsUP (see Definition 1 below) for all $\epsilon$ larger than thresholds that are remarkably close to the $\epsilon_{D}$ above. InAsUP were guessed using trajectory renderings as above.

For large(r) $D$, trajectory renderings are obviously more involved and certainly not so simply useful to detect ergodicity/symmetry breaking. Following standard diagnostics in statistical physics, one can instead use order parameter (OP) empirical estimates. An empirical OP consists of unsigned averages over consecutive iterates of an asymmetry-related observable, which is designed to suggest the existence of asymmetry ergodic components when positive. In order to establish failure of sign-flip symmetry, we use the ”central” coordinate $x_{\lceil\frac{D}{2}\rceil}$ as an observable (Fig. 2); other estimates based on different quantifiers e.g. any single coordinate, two/all coordinate mean values, etc., all yield similar plots with identical bifurcation values (data not shown). Finite-time effects are accounted for by superimposing results from averages over increasing numbers of iterates. Similarly, dependence on the initial condition is evaluated using multiple runs based upon randomly drawn inputs.

In agreement with ergodicity, the OP in Fig. 2 vanishes at small coupling, for every $D$. However, as soon as $\epsilon$ exceeds some $\epsilon_{D}$, this quantity takes on positive values for a positive fraction of initial conditions. This was observed for all investigated values of the dimension, from $D=2$ up to $D=200$. For $D=2$ and $3$, the data are consistent with the phase space plots of Fig. 1; the emergence of positive OP coincides with the appearance of asymmetric ergodic components.444Surprisingly, for $D=3$, the additional asymmetric group emerging at $\epsilon\simeq 0.437>\epsilon_{3}$ in Fig. 1 – signs of the corresponding transition are barely visible on Fig. 2 – shows OP values that are very close, if not identical, to the primary asymmetric group.

In addition to clear evidence of ergodicity breaking, Fig. 2 reveals various interesting phenomenological features of the maps $G_{D,\epsilon}$ dynamics. At first, the bifurcation diagrams show characteristics that are specific to the parity of $D$. For $D$ odd, the bifurcation values $\epsilon_{D}\sim 0.42$ appear to be almost insensitive to $D>3$, and OP estimates are localized around a single ($\epsilon$-dependent) value. Moreover, the fraction of initial conditions that yield non-zero estimates decreases with $D$, making it more difficult to collect marked evidence of ergodicity breaking. For $D$ even, $\epsilon_{D}$ increases with $D$ and seems to approach a limit value $\epsilon_{\ast}<\frac{1}{2}$. OP estimates appear to be uniformly distributed between 0 and the ($\epsilon$-dependent) maximal value, making it difficult to discriminate limit values from finite size fluctuations.

Furthermore, maximal OP values for $\epsilon>\epsilon_{D}$ show a linear $\epsilon$-dependence of tiny slope, which eventually vanishes for large $D$. Also, these maxima decrease as $D$ increases, and asymptotically behave as $\frac{1}{2D}$ for $D$ even (resp. $\frac{1}{D}$ for $D$ odd), see Fig. 3. Accordingly, to unambiguously distinguish asymmetry from short time fluctuations requires longer averages when $D$ increases.555In particular, $T=10^{4}$ averages do not suffice to identify the emergence of positive OP for dimensions $D\geq 70$. This issue, when combined with linear increase of the dimension of the variables in the iterations, substantially increases the computation time required for conclusive evidence. For instance, to obtain the $100\times 100$ averages over $T=10^{6}$ iterates in each of the pictures for $D\geq 70$ required running times of $\sim 100\text{h}$ on a $2.4$ GHz multi-processor computer (compared to $\sim 1\text{min}$ for $D=2$). Therefore, to show ergodicity breaking for dimensions that are commensurate with realistic physical systems appears to be a considerable challenge.

4 Computer proof of ergodicity breaking

The previous section evidences of emergence of asymmetric Lebesgue ergodic components call for rigorous confirmation, in particular to exclude transient effects and other computer round-off shortcomings that may impact numerical simulations and finite-time estimates. A computer-based rigorous proof of ergodicity breaking is presented in this section.

For piecewise affine and expanding maps such as $G_{D,\epsilon}$, any forward invariant finite union of (convex) polytopes must support an absolutely continuous invariant measure [32]. Therefore, in order to prove existence of asymmetric Lebesgue ergodic components, it suffices to obtain such unions of polytopes that are disjoint from their image under $-\text{\rm Id}_{\mathbb{T}^{D}}$.666Throughout this section, a set $S$ is said to be asymmetric iff $S\cap-\text{\rm Id}_{\mathbb{T}^{D}}(S)=\emptyset$. Accordingly, we shall rely on the following notion.

Definition 1**.**

Let $D,\epsilon$ and a finite union of (non-empty) polytopes $P\subset\mathbb{T}^{D}$ be given. Then $P$ is said to be an InAsUP (acronym for Invariant Asymmetric Union of Polytopes) if it satisfies the following conditions

[TABLE]

As already mentioned, analytic proofs of existence of InAsUP have been established for $D=2,3$. However, the proofs rely on observations of trajectories and would be hardy generalisable to large values of $D$. Instead, a fully computational approach is proposed below, which consists of an algorithm designed to generate InAsUP for arbitrary $D$.

4.1 Principles of the algorithm and exact computer results

In few words, the algorithm simply consists in applying repeated iterations of $G_{D,\epsilon}$ to an asymmetric initial polytope.777See subsection 4.2.1 below for more details and a pseudo-code is given in Appendix E). It terminates when the resulting union set becomes $G_{D,\epsilon}$-invariant, or prematurely stops if the set under construction happens to intersect it image under $-\text{\rm Id}_{\mathbb{T}^{D}}$.

As analytic proofs did, the algorithm optimizes its input by using dynamical information from simulations. Initial polytopes are chosen among cylinder sets that are given by symbolic codes associated with empirical trajectories of positive OP (subsection 4.2.2).

That initial polytopes are cylinder sets is crucial when $\epsilon\in\mathbb{Q}$ because all algorithmic calculations then involve rational numbers only (subsection 4.2.4). Using exact computer arithmetics on such numbers888In particular, we use the GNU arithmetic library GMP [15]., this implies that when the construction completes, the resulting asymmetric invariant set must be a genuine InAsUP. In other words, when the construction completes, the computer provides a rigorous proof of ergodicity breaking for the pair $(D,\epsilon)$ under consideration.

To employ exact arithmetic is an important feature of the InAsUP algorithmic construction. Indeed, analytic proofs for $D=2,3$ have revealed that some polytope facets are exactly mapped onto other ones, even when $\epsilon$ is an irrational number. Hence, some analytic cancellations must take place in the construction, which the algorithm would have to carefully monitor, if it dealt instead with floating-point arithmetic. In short terms, in using floating-point arithmetic, control of round-off errors does not suffice to rigorously assert invariance; hence the choice of exact arithmetic here.

Results on exact computer InAsUP constructions are summarized in the following formal statement.

Claim 2**.**

For all values of the pair $(D,\epsilon)$ in Table 1, the map $G_{D,\epsilon}$ has an InAsUP.

To our best knowledge, these results provide the first proof of ergodicity breaking of $G_{D,\epsilon}$ for $D\geq 4$, and hence of the coupled maps $F_{N,\epsilon}$ for $N\geq 5$. For $D=2,3$, the results confirm the previous analytic proof conclusions. For $D\leq 4$, InAsUP have also been obtained for other values of $\epsilon$ (data not shown). As Fig. 2 suggests, we expect InAsUP to exist for every $\epsilon>\epsilon_{D}$.

In addition, Table 1 provides success ratio statistics against the length of the initial cylinder. Clearly, the larger the length is, the higher a success ratio results. This suggests that it suffices to choose a sufficiently small neighborhood of any point of some asymmetric trajectory to generate an InAsUP. Not only ergodicity breaking is not a spurious numerical illusion, but transient behaviors and round-off errors in simulations (Fig. 2) do not impact this phenomenon.

For completeness, Table 1 also provides statistics about InAsUP cardinality and related CPU computation times (viz. typical value and maximum).999As explained in the end of subsection 4.2.2, polytopes in InAsUP may overlap. InAsUP cardinality thus makes no obvious sense other than justifying the measured CPU times. These quantities show substantial variations upon initial cylinder. However, their variation ranges barely vary with the cylinder length. Furthermore, Table 1 reveals rapid growth in these statistics - which accompanies an exponential growth of the atomic partition cardinality - as the dimension $D$ increases. Therefore, even though the construction could in principle operate for any $D$, these exploding features, which imply similar demand of computational resources, may prevent the construction to terminate for large $D$.101010In particular, for $D=6$ ($N_{6}=4683$), no construction has terminated, due to insufficient available RAM and CPU time.

4.2 Details of the algorithm for InAsUP construction and its numerical implementation

A copy of the source files of the InAsUP construction code (in C) is available online [13] (and, again, see Appendix E for a pseudo-code). Most of the procedures in the algorithm may appear to be standard. However, to manipulate and to apply geometric and topological operations on polytopes in arbitrary dimension (such as testing inclusion, testing intersection, computing intersection set, etc), seem not so common. Thanks to the isotropic nature of the $G_{D,\epsilon}$, a convenient, dimension-independent, approach to polytope manipulation and their operations can be developed in our setting. The purpose of this section is to provide insights into these specific procedures and their numerical implementation.

For convenience, we regard $G_{D,\epsilon}$ as a mapping from $S_{D}=(-\tfrac{1}{2},\tfrac{1}{2})^{D}$ into $\mathbb{T}^{D}$ (to be combined with the canonical projection from $\mathbb{T}^{D}$ to $S_{D}$), and we denote by

[TABLE]

the partition of $S_{D}$ into polytopes $A_{a}$ - called atoms - on which the mapping is affine (or, equivalently, on which the vector function $B_{D}$ is constant). Important for future purposes, every atom is a polytope whose facets are included in discontinuity planes (or in parallel planes). Given the expression of $B_{D}$, every discontinuity plane is characterized by the condition $h(\sum_{k=i}^{j}x_{k})\in\tfrac{1}{2}+\mathbb{Z}$ for some $i<j$.

4.2.1 Algorithmic procedure

The InAsUP construction algorithm can be sketched as follows:

• $\bullet$

Choose an initial polytope $P^{0}\subset\mathbb{T}^{D}$ such that $P^{0}\cap-\text{\rm Id}_{\mathbb{T}^{D}}(P^{0})=\emptyset$.

• $\bullet$

Compute subsequent iterates $P^{t+1}=G_{D,\epsilon}(P^{t})$ for $t\in\mathbb{N}$ – which consist of unions of polytopes111111More precisely, $P^{t+1}$ consists of a single polytope as long as ($t$ is such that) $P^{t}\subset A_{a}$ for a single $a$. However, as soon as we have $P^{t}\cap A_{a_{i}}\neq\emptyset$ for several $i$, then

$G_{D,\epsilon}(P^{t})=\bigcup_{i}G_{D,\epsilon}|_{A_{a_{i}}}(P^{t}\cap A_{a_{i}})\ \text{mod}\ 0$

viz. $G_{D,\epsilon}(P^{t})$ becomes the union of several polytopes (mod 0), and therefore so must be all subsequent images $P^{t^{\prime}}$ for $t^{\prime}>t$. – until either $P^{t+1}\subset\bigcup_{k=0}^{t}P^{k}\ \text{mod}\ 0$, or some polytope in $-\text{\rm Id}_{\mathbb{T}^{D}}(P^{t+1})$ intersects $\bigcup_{k=0}^{t+1}P^{k}$.

As explained before, the purpose of the second terminating condition is to interrupt the construction when guarantee fails that the constructed sets will be asymmetric. Otherwise, when the algorithm terminates in the first instance, then both resulting $\bigcup_{k=0}^{t}P^{k}$ and $-\text{\rm Id}_{\mathbb{T}^{D}}(\bigcup_{k=0}^{t}P^{k})$ must be disjoint InAsUPs, as desired.

The map $G_{D,\epsilon}$ does not show any basic property, such as Markov partition [22], that would ensure the algorithm to always terminate in finite time. However, this is the case of all results in Table 1, viz. when the construction does not succeed, intersection with symmetric image is always the cause.

4.2.2 Adapted polytopes and their vector representation

A major aspect of the numerical implementation is to establish a suitable consideration of polytopes and their dynamics. This is the purpose of this subsection.

Since the affine part of $G_{D,\epsilon}$ is a multiple of the identity, the map $G_{D,\epsilon}$ itself preserves polytope facets orientation. A moment of reflection then concludes that the algorithm may deal exclusively with convex polytopes whose facets are aligned with discontinuity planes (because the corresponding set is invariant both under forward and backward dynamics). In particular, this is the case when the initial polytope is a cylinder set, viz. $P^{0}=A_{a_{0}\dots a_{\ell}}$ for some word $\{a_{k}\}_{k=0}^{\ell}$ with letters $a_{k}\in\{1,\dots,N_{D}\}$, where $A_{a_{0}\dots a_{\ell}}$ is defined by [22]

[TABLE]

In practice, candidate words $\{a_{k}\}_{k=0}^{\ell}$ for non-empty cylinders are obtained from empirical trajectories (whose OP is positive), by recording successive labels of visited atoms.

Besides, following standards [7], polytopes in arbitrary dimension can be characterized as intersections of half-spaces, using inequality constraints on coordinates. Half-spaces in our context are delimited by discontinuity planes. Hence polytopes in the algorithm will be characterized using inequality constraints on the sums $\sum_{k=i}^{j}x_{k}$ of coordinates. Accordingly, given a vector

[TABLE]

we define the polytope $P_{m}$ by121212Strict inequalities in this definition are chosen on purpose. Indeed, excluding polytope facets is a convenient way to exclude discontinuities from InAsUP construction.

[TABLE]

(provided that this set is not empty). Overdetermination in this expression is a convenient way to capture in a unique formal expression, all possible types of polytopes that can occur in the construction. In particular, while the number of atom facets may vary, every atom can be expressed as a polytope $P_{m}$ ie. for every $A_{a}\subset S_{D}$, we have $A_{a}=P_{m_{a}}$ for some $m_{a}\in(\frac{1}{2}+\mathbb{Z})^{D(D+1)}$.

Not only polytopes can be represented by vectors (so that the polytope dynamics reduces to that of a $D(D+1)$-dimensional dynamical system), but some of their operations can be expressed in terms of vectors. The following operations are of special interest for our purpose. In all lines below, the indices $i,j$ run over all pairs such that $1\leq i\leq j\leq D$.

• $\bullet$

Inclusion: $P_{m}\subset P_{m^{\prime}}$ if $\underline{m^{\prime}}_{i,j}\leq\underline{m}_{i,j}$ and $\overline{m}_{i,j}\leq\overline{m^{\prime}}_{i,j}$.

• $\bullet$

Intersection: $P_{m}\cap P_{m^{\prime}}=P_{m\cap m^{\prime}}$ where $\underline{(m\cap m^{\prime})}_{i,j}=\max\{\underline{m}_{i,j},\underline{m^{\prime}}_{i,j}\}$ and $\overline{(m\cap m^{\prime})}_{i,j}=\min\{\overline{m}_{i,j},\overline{m^{\prime}}_{i,j}\}$.

• $\bullet$

Dilation: For every $a>0$, we have $a\text{Id}_{\mathbb{T}^{D}}(P_{m})=P_{m^{\prime}}$, where $\overline{\underline{m^{\prime}}}_{i,j}=a\overline{\underline{m}}_{i,j}$.

• $\bullet$

Translation: For every $x\in\mathbb{R}^{D}$, we have $P_{m}+x=P_{m^{\prime}}$ where $\overline{\underline{m^{\prime}}}_{i,j}=\overline{\underline{m}}_{i,j}+\sum_{k=i}^{j}x_{k}$.

• $\bullet$

Symmetry: $-P_{m}=P_{\Sigma(m)}$ where $\underline{\Sigma(m)}_{i,j}=-\overline{m}_{i,j}$ and $\overline{\Sigma(m)}_{i,j}=-\underline{m}_{i,j}$.

The union operation is however not so convenient. In particular, depending on $m$ and $m^{\prime}$, there might not exist $m^{\prime\prime}\in\mathbb{R}^{D(D+1)}$ such that

[TABLE]

Therefore, when testing that $P^{t+1}\subset\bigcup_{k=0}^{t}P^{k}\ \text{mod}\ 0$, the algorithm can only test if for every $P_{m}\subset P^{t+1}$, there exists $P_{m^{\prime}}$ in the collection $\bigcup_{k=0}^{t}P^{k}$ such that $P_{m}\subset P_{m^{\prime}}$.

While the intersection property ensures immediate detection of asymmetry failure, a consequence of this constrained inclusion test is that the algorithm may continue to run while InAsIUP construction has been completed. In particular, polytopes in the final union may overlap and the corresponding cardinality might be unnecessarily excessive.

Polytopes overlap can however be reduced by chopping off pieces. Indeed one can test if a polytope $P_{m}\subset P^{t+1}$ writes $P_{m}=\bigcup_{\ell}P_{m_{\ell}}\ \text{mod}\ 0$ with $P_{m_{\ell}}\subset P_{m^{\prime}}\subset\bigcup_{k=0}^{t}P^{k}$ for some $\ell$, so that the construction can only retain the complementary $P_{m_{\ell}}$ that do not belong to any $P_{m^{\prime}}$ in the existing union. Retaining smaller polytopes makes it more likely that the inclusion $P^{t+1}\subset\bigcup_{k=0}^{t}P^{k}\ \text{mod}\ 0$ holds for $t$ smaller, and thus that the final union cardinality is smaller. Such a diminution has been observed in practice, accompanied with substantial reduction of computer resources when $D>3$.

For simplicity, we opted to test if $P_{m}\subset P^{t+1}$ decomposes as the union of two pieces, namely

[TABLE]

with $P_{m_{\text{in}}}=P_{m}\cap P_{m^{\prime}}$ for some $P_{m^{\prime}}\subset\bigcup_{k=0}^{t}P^{k}$ (so that it suffices to retain $P_{m_{\text{out}}}$), see Appendix D for mathematical foundations and criteria.131313Actually, a more elaborated version of the algorithm also tests the three-piece decomposition

$P_{m}=P_{m_{\text{in}}}\cup P_{m_{\text{out},1}}\cup P_{m_{\text{out},2}}\ \text{mod}\ 0$

in the simplest case where $P_{m}$ intersects $P_{m^{\prime}}$ transversally, so that the $P_{m_{\text{out},i}}$ result from chopping off along parallel faces.

4.2.3 Polytope constraint optimization

While expression (4.1) is convenient, it does not implies that $P_{m}$ is not empty, nor uniqueness of the constraining vector $m$. Indeed, $P_{m}$ can be unambiguously specified and yet, as many as $D(D-1)$ constraints may remain inactive (Fig. 4). Inactive constraints [7] are problematic in an automated numerical implementation because they may yield spurious polytopes.

These issues call for constraint optimization. For the sake of notations, we present an optimization procedure in a slightly more general setting. The procedure relies on standard Linear Programming arguments.

Let $D\leq E$ be arbitrary integers. Suppose that a collection $\alpha=\left((\alpha_{ij})_{j=1}^{D}\right)_{i=1}^{E}$ of distinct vectors with non-negative elements is given.141414That is to say, $\alpha_{ij}\geq 0$ and $(\alpha_{ij})_{j=1}^{D}\neq(\alpha_{i^{\prime}j})_{j=1}^{D}$ when $i\neq i^{\prime}$. For any constraint vector $m=(\underline{m}_{i},\overline{m}_{i})_{i=1}^{E}$, consider the polytope $P_{m}^{\alpha}$ defined by

[TABLE]

For each $i\in\{1,\dots,E\}$, consider the set $\Lambda_{i}$ of vectors $(\lambda_{k})_{k=1}^{E}$ with at least $E-D$ vanishing coordinates, which uniquely solve the system of equations

[TABLE]

More precisely, given any $s\in\{1,\dots,D\}$ and $S\subset\{1,\dots,E\}$ of cardinality $s$, if it exists, let $\lambda^{S}$ be the unique solution of the equations obtained from (4.2) by letting $\lambda_{k}=0$ for $k\in\{1,\dots,E\}\setminus S$. The set $\Lambda_{i}$ is made of all such solutions $\lambda^{S}$ when $S$ and $s$ vary.151515In particular, $\Lambda_{i}$ is a finite set which can be obtained by listing all possible sets $S$ and computing, when they exist, the unique solutions $\lambda^{S}$ of the corresponding systems. Moreover, $\Lambda_{i}$ cannot be empty because it contains the canonical vector $(\lambda_{k})_{k=1}^{E}=(\delta_{k,i})_{k=1}^{E}$, where $\delta_{k,i}$ is the Kronecker symbol.

Independently, given $\lambda\in\mathbb{R}$ and a constraint vector $m$, consider the vectors $(\underline{e}_{k}(\lambda,m))_{k=1}^{E}$ and $(\overline{e}_{k}(\lambda,m))_{k=1}^{E}$ respectively defined by

[TABLE]

Optimized constraint vectors, together with existence condition, are given in the following statement.

Lemma 3**.**

Given any constraint vector $m=(\underline{m}_{i},\overline{m}_{i})_{i=1}^{E}$, consider $O(m)=(\underline{O(m)}_{i},\overline{O(m)}_{i})_{i=1}^{E}$ defined by

[TABLE]

Then

• (i)

$P_{m}^{\alpha}$* is not empty iff $\underline{O(m)}_{i}<\overline{O(m)}_{i}$ for all $1\leq i\leq E$.*

• (ii)

If $P_{m}^{\alpha}$ is not empty, then $P_{O(m)}^{\alpha}=P_{m}^{\alpha}$ and all constraints in the definition of $P_{O(m)}^{\alpha}$ are active.

• (iii)

The plane $\sum_{j=1}^{D}\alpha_{ij}x_{j}=\underline{O(m)}_{i}$ (resp. $\sum_{j=1}^{D}\alpha_{ij}x_{j}=\overline{O(m)}_{i}$) defines a face of $P_{m}^{\alpha}$ iff

[TABLE]

The proof, given in Appendix C, is inspired from the Karush-Kuhn-Tucker (KKT) approach to linear programming, an extension of the method of Lagrange multipliers to the case of inequalities constraints [7].

The algorithm extensively relies on Lemma 3 and in particular, replaces $m$ by $O(m)$ every time the intersection of two polytopes is tested or computed.

For the collection $\alpha$ involved in the definition (4.1), the cardinality of the sets $\Lambda_{i}$ does not depend on $i$ but it exponentially increases with $D$, see Table 2. According to the GNU performance analysis tool gprof, the optimization procedure $m\leftarrow O(m)$ is the most consuming resource task of the overall execution, in terms of CPU time.

In order to speed up the process, the results in Table 1 have actually been obtained using the semi-optimization procedure $O^{\prime}$ which maximises/minimises over the subsets $\Lambda^{\prime}_{i}\subset\Lambda_{i}$ of vectors having at most two non-vanishing coordinates. Naturally, the resulting vectors $O^{\prime}(m)$ are no longer optimal (hence some of the $P_{O^{\prime}(m)}$ may be spurious). Therefore, the cardinality of the constructing InAsIUP may be larger than when using $O$ only. However, this deficiency does not seem to affect InAsIUP successful completion in practice, especially because the chopping procedure of Appendix D still applies. More importantly, computation times appear to be reduced by a factor $\sim\tfrac{\#\Lambda_{i}}{\#\Lambda^{\prime}_{i}}$ which, since $\#\Lambda^{\prime}_{i}=D$, is a substantial gain for $D\geq 4$.

4.2.4 Polytope related vector dynamics

Let $B_{a}=B_{D}|_{A_{a}}$ be the expression in atom $A_{a}$ of the vector involved in the constant part of $G_{D,\epsilon}$. In each atom, $G_{D,\epsilon}$ regarded as acting on polytopes $P_{m}$ induces a mapping $\Gamma_{D,\epsilon,a}$ on vectors $m$, via the relation $G_{D,\epsilon}|_{A_{a}}(P_{m})=P_{\Gamma_{D,\epsilon,a}(m)}$. Explicitly, we have

[TABLE]

Moreover, the atomic decomposition

[TABLE]

where $P_{m\cap m_{a}}=P_{m}\cap P_{m_{a}}$ and $P_{m_{a}}=A_{a}$, induces an extension of the vector map $\Gamma_{D,\epsilon}$ to arbitrary polytopes $P_{m}\subset S_{D}$, ie. we have

[TABLE]

where the label index $a$ runs over all labels for which $P_{m\cap m_{a}}$ is not empty. Furthermore, in relation with the previous section, notice that the vector dynamics preserve the optimization procedure, ie. we have

[TABLE]

Notice finally that, since every $B_{a}\in\mathbb{Q}^{D}$, when $\epsilon\in\mathbb{Q}$ a rational number, every map $\Gamma_{D,\epsilon,a}$ has the properties

[TABLE]

viz. the polytope related vector dynamics effectively implies only rational numbers when $\epsilon\in Q$.

5 Discussion

The literature on the deterministic dynamics of collective systems mostly describes phenomenological changes related to pattern formation or transition to synchrony. These changes typically correspond to bifurcations of steady states or periodic solutions [1, 11, 10, 28]. These transitions – which may also imply a reduction in the phase space to lower dimensional subspaces [2] – can be regarded as the breaking of ergodicity of an atomic or singular measure. Instead, analogues of phase transitions with spontaneous symmetry breaking should involve ACIM, in order to preserve all degrees of freedom of the corresponding chaotic attractor. Rigorously articulated examples of this are few, especially when appeal to a underlying ad hoc random process possessing the desired property is excluded.

In this paper, we have first provided (complementary) numerical evidence for symmetry breaking of a chaotic attractor of full dimension, in a model of a collective system of interacting individuals (whose Markov partition remains elusive). Furthermore, we have developed and benchmarked an algorithm for exact computer proof of the corresponding breaking of ACIM ergodicity. Even though the algorithm could only terminate for systems with limited number of individuals (due to limitations of computational resources), it indicates that phase transitions can be rigorously proven in a purely deterministic setting, without any reference to statistical mechanics.

In order to improve their physical relevance, such features should be confirmed for systems with larger numbers of individuals and for more realistic models. For our systems, such confirmation might require algorithmic improvements. For instance, one may rely on parallel implementation – even though the iterative construction of an invariant set is an intrinsically sequential process – on automated decomposition into unions of non-overlapping polytopes, and on the use of optimized libraries for rational arithmetic, as well as dedicated computational resources.

Besides the existence of InAsUP, alternative proofs of the breaking of ACIM ergodicity could be obtained using spectral properties of the transfer operator associated with the coupled map system. Since this operator governs the dynamics of the densities associated with measures, it suffices to prove that it acquires multiple fixed points [18] as the coupling intensity increases. Likewise, one could aim at delineating (coupling-dependent) Markov partitions, which are compatible with ergodicity at small coupling and simultaneously imply splitting into asymmetric ergodic components as interactions become strong. To our best knowledge, these approaches have not been considered in the literature.

In any case, while transitions in statistical mechanics only occur in the thermodynamic limit, $N\to+\infty$ – especially because irreducible Markov chains, which govern the dynamics for finite $N$, must be uniquely ergodic – this paper, together with [12, 30, 31], shows that without Markov chain considerations, deterministic models do not require to consider this limit, and can display similar non-ergodic behaviors and symmetry breaking in finite dimension. This feature is particularly interesting for the modelling of real particle systems, which must be finite-dimensional.

Acknowledgments

I am grateful to P. Bálint, J. Buzzi, V. Perchet, F. Sélley and L-S. Young for scientific discussions, to Y. Legrandgérard, L. Ollivier and D. Simon for computational insights and to P. Bálint, N. Cuneo and G. Francfort for a critical reading of the manuscript and thoughtful suggestions of improvements.

Appendix A Symmetry breaking for $G_{2,\epsilon}$

The map $G_{2,\epsilon}$’s explicitly expression is given by $G_{2,\epsilon}=2(1-\epsilon)\text{\rm Id}_{\mathbb{T}^{2}}+\tfrac{2\epsilon}{3}B_{2}\ \text{mod}\ 1$, where

[TABLE]

$G_{2,\epsilon}$ commutes with each transformation in the natural representation of $\mathbb{Z}_{2}\times S_{3}$ on $\mathbb{T}^{2}$. This group can be generated by the sign flip $-\text{Id}_{\mathbb{T}^{2}}$ and the transformations $\sigma_{213},\ \sigma_{132}$ and $\sigma_{321}$, which are induced by the transpositions in $S_{3}$ (NB: subscripts here denote transposition images of the ordered list $\{1,2,3\}$ - abbreviated as 123). The sign flip writes

[TABLE]

and it corresponds to the central symmetry wrt $(\frac{1}{2},\frac{1}{2})$ in the square $(0,1)^{2}$. Similarly, we have

[TABLE]

which corresponds to the orthogonal reflection wrt to the anti-diagonal $x_{1}+x_{2}=1$. Moreover, $\sigma_{213}$ (resp. $\sigma_{132}$) is defined by

[TABLE]

is the reflection wrt the axis $x_{1}=0$ along the direction $x_{1}+2x_{2}=\text{cst}$ (resp. wrt the axis $x_{2}=0$ along the direction $2x_{1}+x_{2}=\text{cst}$).

The representation on $\mathbb{T}^{2}$ of $\mathbb{Z}_{2}\times S_{3}$ itself can be listed as follows

[TABLE]

Clearly, all ergodic components for $\epsilon>\epsilon_{2}$ break the sign-flip symmetry (Fig. 1). More precisely, each component remains invariant under the action of one transformation above and is mapped onto another component under any other transformation. In particular, the chartreuse and dark green components are invariant under $\sigma_{213}$, Blue and light green components are invariant under $\sigma_{132}$ and red and orange components are invariant under $\sigma_{321}$.

Appendix B Symmetry breaking for $G_{3,\epsilon}$

The map $G_{3,\epsilon}$ writes $G_{3,\epsilon}=2(1-\epsilon)\text{\rm Id}_{\mathbb{T}^{3}}+\tfrac{\epsilon}{2}B_{3}\ \text{mod}\ 1$, where

[TABLE]

As before, any transformation in the representation of $S_{4}$ on $\mathbb{T}^{3}$ can be obtained as a composition of 2-element-permutation representations, whose expressions are given by (NB: the $(\text{mod}\ 1)$ are not specified for the sake of space)

[TABLE]

The ergodic components that emerge at $\epsilon=\epsilon_{3}$ all break the sign-flip symmetry $-\text{Id}$ in $\mathbb{T}^{3}$. As before, these 6 components are only partly asymmetric; they remain invariant under the action of 7 transformations in $\mathbb{Z}_{2}\times S_{4}$, and they are exchanged under other transformations. In particular, the blue component is invariant under

[TABLE]

and, obviously, under the composition $\sigma_{3214}\circ\sigma_{1432}=\sigma_{1432}\circ\sigma_{3214}$.

As for the 8 components emerging at $\epsilon\simeq 0.437$, their asymmetries are stronger than for the previous components, as they remain invariant under only 5 transformations. For instance, the fuschia component is invariant under

[TABLE]

and, obviously, under the compositions $\sigma_{4231}\circ\sigma_{2134}=\sigma_{2134}\circ\sigma_{4231}=\sigma_{2134}\circ\sigma_{1432}=\sigma_{1432}\circ\sigma_{4231}$ and $\sigma_{1432}\circ\sigma_{2134}=\sigma_{4231}\circ\sigma_{1432}$.

Appendix C Proof of Lemma 3

We first establish the expression of $\overline{O(m)}_{i}$. A similar analysis yields the expression of $\underline{O(m)}_{i}$. Let $i\in\{1,\dots,E\}$ and suppose that we aim to maximize the quantity $\sum_{j=1}^{D}\alpha_{ij}x_{j}$ given the $2E$ linear inequality constraints

[TABLE]

This problem exactly fits the KKT approach to linear programming [7]. The corresponding KKT conditions state in particular that every maximizer $(x_{j}^{\ast})_{j=1}^{D}$ must cancel the gradient of the following Lagrangian

[TABLE]

for a unique pair of vectors $(\gamma_{k})_{k=1}^{E},(\gamma^{\prime}_{k})_{k=1}^{E}$ with non-negative components. In other words, this pair of vectors must satisfy the equation

[TABLE]

The complementary slackness conditions in the KKT setting then state that we must also have

[TABLE]

Now, the conditions $\underline{m}_{k}<\overline{m}_{k}$ imply that the constraints $\sum_{j=1}^{D}\alpha_{kj}x_{j}^{\ast}\leq\overline{m}_{k}$ and $-\sum_{j=1}^{D}\alpha_{kj}x_{j}^{\ast}\leq-\underline{m}_{k}$ cannot be simultaneously active, viz. we must have $\gamma_{k}\gamma^{\prime}_{k}=0$ for all $k$. A single multiplier $\lambda_{k}=\gamma_{k}-\gamma^{\prime}_{k}$ results for each $k$ and equation (C.2) implies that the vector $(\lambda_{k})_{k=1}^{E}$ must solve equation (4.2).

Equation (C.3) also implies that for $\lambda_{k}\neq 0$ we have

[TABLE]

Together with equation (4.2), this equality yields

[TABLE]

which is exactly the expression involved in the definition of $\overline{O(m)}_{i}$.

When $E>D$, the system (4.2) is underdetermined. Hence it may have infinitely many solutions. The following considerations show that we may only retain solutions in $\Lambda_{i}$.

Consider the partition of $\mathbb{R}^{E}$ into orthants in the interior of which the sign of every coordinate $\lambda_{k}$ is constant. Choose any orthant that contains a family of (infinitely many) solutions to (4.2) in its interior.161616If no such orthant exists, then proceed with similar considerations in orthant boundaries, as indicated in the next paragraph. Those coordinates $\lambda_{k}$ with identical sign ($E>2$) must have opposite variations when varying the solution in the family. However, the functional $(\lambda_{k})_{k=1}^{E}\mapsto\sum_{k=1}^{E}\lambda_{k}\overline{e}_{k}(\lambda_{k},m)$ is linear (and hence with constant gradient) inside the orthant. Therefore, it certainly reaches its maximum on the boundary of the orthant, ie. when at least one coordinate of $(\lambda_{k})_{k=1}^{E}$ vanishes.

Repeating the reasoning inside orthant boundaries, which are identified with $\mathbb{R}^{E-1}$, and then for all subsequent subspaces $\mathbb{R}^{E^{\prime}}$ for decreasing $E^{\prime}$, from $E-2$ to $D+1$, viz. as long as the resulting systems remain underdetermined, we conclude that we may only consider solutions of (4.2) with $E-D$ vanishing coordinates.

Yet, the system with $D$ unknowns/equations may also be degenerate. In this case, one can remove one superfluous equation from (4.2) and repeat the previous argument to conclude that the maximum must occur for a vector with $E-D+1$ vanishing coordinates. Repeating this reasoning as long as a degeneracy occurs, we may eventually get a system of 2 unknowns/equations. This system may be problematic if degenerate and the two coordinates have opposite signs. In fact, the only problematic scenario is if the functional $\sum_{k=1}^{E}\lambda_{k}\underline{e}_{k}(\lambda_{k},m)$ increases when the positive coordinates increases (otherwise the maximum is reached when one coordinate vanishes). In this case, the maximum would be $+\infty$, which is certainly impossible. This analysis concludes that $\overline{O(m)}_{i}$ maximizes the quantity $\sum_{j=1}^{D}\alpha_{ij}x_{j}$, under the constraints (C.1).

As a consequence, if $P_{m}^{\alpha}\neq\emptyset$, then for every $x\in\overline{P_{m}^{\alpha}}$ (the closure of $P_{m}^{\alpha}$), we must have

[TABLE]

viz. $\overline{P_{m}^{\alpha}}\subset\overline{P_{O(m)}^{\alpha}}$ (and then $P_{O(m)}^{\alpha}\neq\emptyset$). On the other hand, that the canonical vector $(\delta_{k,i})_{k=1}^{E}$ belongs to $\Lambda_{i}$ implies

[TABLE]

hence $\overline{P_{O(m)}^{\alpha}}\subset\overline{P_{m}^{\alpha}}$. It results that $P_{O(m)}^{\alpha}=P_{m}^{\alpha}$ when $P_{m}^{\alpha}\neq\emptyset$.

Now, the complementary slackness conditions (C.3) imply that the maximum $\overline{O(m)}_{i}$ and minimum $\underline{O(m)}_{i}$ are given by combinations of values at active constraints. Besides, we must have $\overline{O(m)}_{i}=\overline{m}_{i}$ and $\underline{O(m)}_{i}=\underline{m}_{i}$ when the constraints at $i$ are active, by the definitions of $\overline{O(m)}_{i}$ and $\underline{O(m)}_{i}$. It follows that the values of $\overline{O(m)}_{i}$ and $\underline{O(m)}_{i}$ do not change when we replace $m$ by $O(m)$ in the constraints (C.1) (viz. $O$ is a projection operator) and the complementary slackness conditions imply that all constraints must be active in the updated optimisation problem, ie. all constraints in the definition of $P_{O(m)}^{\alpha}$ must be active.

Moreover, that $\overline{O(m)}_{i}$ and $\underline{O(m)}_{i}$ are coordinate extremal values in $P_{O(m)}^{\alpha}=P_{m}^{\alpha}$ immediately imply that

[TABLE]

is an existence condition for $P_{m}^{\alpha}$.

Proof of item (iii). If a maximizer $(x_{j}^{\ast})_{j=1}^{D}$ satisfies

[TABLE]

then, by continuity, the inequality

[TABLE]

holds for all points in the intersection of the plane $\sum_{j=1}^{D}\alpha_{ij}x_{j}=\overline{O(m)}_{i}$ with a sufficiently small neighborhood of $(x_{j}^{\ast})_{j=1}^{D}$; hence this plane defines a facet of $P_{m}^{\alpha}$.

On the opposite, if

[TABLE]

then let $\delta=(\delta_{i})_{i=1}^{D}$ be a perturbation of the maximiser $x^{\ast}$ in the plane $\sum_{j=1}^{D}\alpha_{ij}x_{j}=\overline{O(m)}_{i}$, ie. such that $\sum_{j=1}^{D}\alpha_{ij}\delta_{j}=0$. Then for those $(\lambda_{k})$ for which the previous equality holds, we must also have

[TABLE]

Since the vectors $(\alpha_{kj})_{j=1^{D}}$ are all distinct, this means that $\delta$ must be limited to a set of co-dimension $1+\#\{k\ :\ \lambda_{k}\neq 0\}>1$. Therefore, it can certainly not span a facet of $P_{m}^{\alpha}$. The proof of Lemma C is complete. $\Box$

Appendix D Mathematical statements related to polytope chopping tests

The setting in this section assumes that two polytopes $P_{m}^{\alpha}$ and $P_{m^{\prime}}^{\alpha}$ (where the constraint vectors $m$ and $m^{\prime}$ are assumed to be optimized) intersect but are not contained in one another. We aim to determine (simple) conditions so that $P_{m}^{\alpha}\setminus P_{m\cap m^{\prime}}^{\alpha}$ consists of either one or two polytopes defined by inequality constraints of the same type. We shall rely on the following set of indices

[TABLE]

Claim 4**.**

In the current setting, assume that171717We have $i^{-},i^{+}\in S$ and $i^{-}$ (resp. $i^{+}$) exists only when $\underline{m^{\prime}}_{i^{-}}$ (resp. $\overline{m^{\prime}}_{i^{+}}$) defines a facet of $P_{m^{\prime}}^{\alpha}$.

• •

a unique $i^{-}$ exists such that $\underline{m}_{i^{-}}<\underline{m^{\prime}}_{i^{-}}$ and $\underline{O(m\cap m^{\prime})}_{i^{-}}=\underline{m^{\prime}}_{i^{-}}$

• •

and/or, a unique $i^{+}$ exists such that $\overline{m^{\prime}}_{i^{+}}<\overline{m}_{i^{+}}$ and $\overline{O(m\cap m^{\prime})}_{i^{+}}=\overline{m^{\prime}}_{i^{+}}$.

Then, the intersection $P_{m\cap m^{\prime}}^{\alpha}:=P_{m}^{\alpha}\cap P_{m^{\prime}}^{\alpha}$ can be characterized as follows

[TABLE]

The numerical algorithm relies on the following consequence to define complementary piece(s) when an index of type $i^{-}$ or $i^{+}$ appears in testing intersections.

Corollary 5**.**

If $i^{-}$ exists, let the constraint vector $m^{-}$ be defined by

[TABLE]

and similarly, let $m^{+}$ be defined by

[TABLE]

when assuming $i^{+}$ exists. We have

[TABLE]

where $P_{m^{-}}^{\alpha}$ (resp. $P_{m^{+}}^{\alpha}$) has to be replaced by $\emptyset$ if $i^{-}$ (resp. $i^{+}$) does not exist.

The proof of the Corollary is immediate. The definitions of $m^{-}$ and $m^{+}$ immediately follow from the characterizations of $m\cap m^{\prime}$ and $i^{-},i^{+}$. That the corresponding polytopes are non-empty is a consequence of the fact that $m$ and $m^{\prime}$ are optimized (so that there exists $x^{-}\in P_{m^{-}}^{\alpha}$ such that $x_{i^{-}}^{-}=\underline{m}_{i^{-}}$ and $x^{+}\in P_{m^{+}}^{\alpha}$ such that $x_{i^{+}}^{+}=\underline{m}_{i^{+}}$).

The Claim remains valid if one considers the semi-optimization procedure $O^{\prime}$ (obtained by replacing the $\Lambda_{i}$ by $\Lambda^{\prime}_{i}$). However $P_{m^{-}}^{\alpha}\neq\emptyset$ and $P_{m^{+}}^{\alpha}\neq\emptyset$ can no longer be granted. So in using $O^{\prime}$ instead of $O$ when testing chopping, the algorithm may include spurious polytopes in the constructing collection.

Proof of the Claim. Consider all possible cases of ordering of the coordinates of $m$ and $m^{\prime}$ under the assumption $m\cap m^{\prime}\neq\emptyset$. Then, from the definition of $m\cap m^{\prime}$ in Section 4.2.2, we must have

[TABLE]

Assume now that $i\in S$, otherwise there is nothing to prove. Consider first the case $\underline{m}_{i}<\underline{m^{\prime}}_{i}<\overline{m}_{i}\leq\overline{m^{\prime}}_{i}$. Then we must have $i\neq i^{+}$ and we can let $\overline{(m\cap m^{\prime})}_{i}=\overline{m}_{i}$. So,

• •

either $i=i^{-}$ and we get $\underline{(m\cap m^{\prime})}_{i^{-}}=\underline{m^{\prime}}_{i^{-}}$ from the definition of $i^{-}$.

• •

or $i\neq i^{-}$ and then we must have $\underline{m}_{i}<\underline{m^{\prime}}_{i}<\underline{O(m\cap m^{\prime})}_{i}$. Obviously, $\underline{m^{\prime}}_{i}$ cannot define an active constraint of $P_{m\cap m^{\prime}}^{\alpha}$ and we may set $\underline{(m\cap m^{\prime})}_{i}=\underline{m}_{i}$ without affecting this set.

Now, the second case $\underline{m^{\prime}}_{i}\leq\underline{m}_{i}<\overline{m^{\prime}}_{i}<\overline{m}_{i}$ can be treated similarly. The third case is even simpler because either $i=i^{-}=i^{+}$ and then we get $(\underline{m^{\prime}}_{i},\overline{m^{\prime}}_{i})$, or $i\in\{i^{-},i^{+}\}$ and $i^{+}\neq i^{-}$ and we get one of the previous cases. $\Box$

Appendix E Algorithm for InAsUP construction

Declarations, initialisations

  Choose $D,\epsilon$

  Import vectors $(\lambda_{k})_{k=1}^{\frac{D(D+1)}{2}}\in\bigcup_{i=1}^{D}\Lambda_{i}$ from output of automated Mathematica code.

  Import symbolic word $a_{0}\cdots a_{\ell}$ of initial cylinder, from symbolic sequence of empirical trajectory simulated using Mathematica code.

Preliminaries: Computations of atom constraint vectors $m_{a}$ and constant part of $\Gamma_{D,\epsilon,a}$

  Set $a\leftarrow 0$

  List all possible vectors $\left\{h\left(\sum_{k=i}^{j}x_{k}\right)\right\}_{1\leq i\leq j\leq D}$ for $x\in S_{D}$.

  for each vector do

     $m\leftarrow\left((\underline{m}_{i,j},\overline{m}_{i,j})_{i=1}^{j}\right)_{j=1}^{D}$ where $\underline{m}_{i,j}=h\left(\sum_{k=i}^{j}x_{k}\right)-\frac{1}{2}$ and $\overline{m}_{i,j}=h\left(\sum_{k=i}^{j}x_{k}\right)+\frac{1}{2}$

     Apply optimization, ie. $m\leftarrow O(m)$

     Store constraint vector, ie. $m_{a}\leftarrow m$

     Store constant part vector, ie. $B_{a}\leftarrow\frac{2\epsilon}{D+1}\left\{\sum_{k=i}^{j}B_{a,k}\right\}_{1\leq i\leq j\leq D}$

     Set $a\leftarrow a+1$

  end for

Implementation

  Replace $O$ by the semi-optimized procedure $O^{\prime}$, to be used in all future computations.181818For the sake of brevity, optimizations $m\leftarrow O^{\prime}(m)$ are not indicated in the rest of the algorithm. However, they are applied every time the intersection of two polytopes is tested or computed and when chopping into pieces is tested or computed.

  Compute constraint vector $m_{a_{0}\cdots a_{\ell}}$ of initial cylinder recursively, ie. $m_{a_{0}\cdots a_{\ell}}=m_{a_{0}}\cap\Gamma_{D,\epsilon,a_{0}}^{-1}(m_{a_{1}\cdots a_{\ell}})$.

  Set $\text{UP}\leftarrow m_{a_{0}\cdots a_{\ell}}$ and $\text{NP}\leftarrow m_{a_{0}\cdots a_{\ell}}$.

  if $m_{a_{0}\cdots a_{\ell}}\cap\Sigma(m_{a_{0}\cdots a_{\ell}})\neq\emptyset$ then

     Exit implementation (Asymmetry fails)

  end if

  while $\text{NP}\neq\emptyset$ $\&$ asymmetry holds $\&$ $\#\text{UP}<\text{MaxCard}$ (Main Loop) do

     for each $m\in\text{NP}$ do

        Compute image $m^{\prime}=\Gamma_{D,\epsilon,a}(m)$.

        Apply projection $\mathbb{T}^{D}\to S_{D}$. May chop $m^{\prime}$ into several pieces (when $m^{\prime}$ spans across boundaries of $S_{D}$).

        for each piece of $m^{\prime}$ and each symbol $a$ do

           if $m^{\prime}\cap m_{a}\neq\emptyset$ then

              Set $m^{\prime}\leftarrow m^{\prime}\cap m_{a}$

              if $m^{\prime}\cap\Sigma(\text{UP})\neq\emptyset$ $\parallel$ $m^{\prime}\cap\Sigma(\text{EP})\neq\emptyset$ then

                 Exit implementation (Asymmetry fails)

              else

                 for each $m^{\prime\prime}\in(\text{UP}\cup\text{EP})\cap m_{a}$191919In practice, this loop consists of two distinct loops; one for $m^{\prime\prime}\in\text{UP}\cap m_{a}$ followed by one for $m^{\prime\prime}\in\text{EP}\cap m_{a}$. do

                    if $m^{\prime}\cap m^{\prime\prime}\neq\emptyset$ but $m^{\prime}\not\subset m^{\prime\prime}$ then

                       if $m^{\prime}=m^{-}\cup m^{\prime\prime}\cup m^{+}$ ($m^{-}\neq\emptyset$ $\parallel$ $m^{+}\neq\emptyset$) then

                          $m^{\prime}\leftarrow(m^{-}\cup m^{+})$ (Chopping operation).

                       end if

                    end if

                 end for

                 Set $\text{EP}\leftarrow\text{EP}\cup m^{\prime}$.202020Due to chopping, $m^{\prime}$ may consists of several elements.

              end if

           end if

        end for

     end for

     Set $\text{NP}\leftarrow\text{EP}$, $\text{UP}\leftarrow\text{UP}\cup\text{NP}$ and $\text{EP}\leftarrow\emptyset$.

  end while

  print  $\#\text{UP}$ and CPU time to complete implementation.

Annex: Definitions of functions

  Optimization function $m\leftarrow O(m)$ (and semi-optimization function $m\leftarrow O^{\prime}(m)$).

  Mapping $\Gamma_{D,\epsilon,a}$ in each atom.

  Pre-image $\Gamma_{D,\epsilon,a}^{-1}$ in each atom.

  Compute optimized intersection $O(m\cap m^{\prime})$

  Test $m\cap m^{\prime}\neq 0$

  Compute symmetric image $\Sigma(m)$.

  Test chopping $m$ by $m^{\prime}$ into one or two polytope(s).

  Compute chopped pieces $m=m^{-}\cup(m^{\prime}\cap m)\cup m^{+}$