Lasers, stability, and numbers

Can insight familiar from quantum physics be used to uncover unsuspected features of classical dynamical systems? Is it possible to discover signatures of familiar quantum concepts in the realm of classical dynamics? For the wide class of systems governed by polynomial equations of motion, this paper argues that the concept of entanglement has a purely classical analogon, illustrating the analogy analytically for two paradigmatic systems, in one and in two-dimensions. The analogy emerges from properties of simple equations of motion of lasers, by demonstrating that complete sets of periodic modes may be encoded into a single orbital equation, a mode carrier, which parameterizes the entire set. The open question of mode isomorphism is also briefly addressed in the context above. Carriers provide a fresh technique for tracking and controlling classical dynamical systems, putting the emphasis on the study of orbital equations instead of orbital points.


Introduction
The stability and dynamical behavior of lasers is known to be well-described by rate equations [1][2][3][4][5]. For instance, the top row of figure 1 illustrates stability diagrams computed for a specific section of the control parameter space of a CO 2 with opto-electronic feedback governed by a state-of-the-art sixdimensional set of rate equations [3]. Similarly, figure 1(c) presents a stability phase diagram obtained for a popular three rate equations model describing oscillatory modes of a semiconductor lasers with optical injection [3]. Figure 1(a) displays a standard Lyapunov phase diagram showing the dichotomic separation of periodic modes, represented by dark shadings (negative Lyapunov exponent), and the non-periodic (chaotic) modes, represented by colors (positive exponents). Complementing the Lyapunov diagram, figure 1(b) shows a so-called isospike diagram, where colors record the number of spikes (local maxima) of the laser intensity detected per period of oscillation. Black indicates lack of periodic oscillations, i.e.chaos. For a recent survey on the use of rate equations to describe stability of laser modes, periodic or not, see [3]. Figure 1(d) shows a familiar bifurcation diagram illustrating laser intensity recorded along the dashed white line marked in figure 1(c). Figure 1 displays conspicuous features abundantly observed in other lasers and in paradigmatic oscillators in several disciplines [6]. First, figure 1 shows stability diagrams for undriven (i.e.autonomous) laser systems. This means that their self-similarities result from an intrinsic auto-organization of the laser modes, not from external modulations. Second, the Lyapunov diagrams on the left column reveal clear accumulations of self-similar stability regions of periodic laser oscillations. The number of spikes per period grows without bound as the islands accumulate more and more towards the boundary on the upper right corner. Third, the isospike diagram reveals that the inner boundaries of the 'shrimp-shaped' stability regions [7] are also self-similar and organized around infinite networks of exceptional points [7][8][9] which define the orientation of the islands in control parameter space.
A major obstacle when one tries to address analytically the regularities in figure 1 is the total lack of a theoretical framework to predict or explain the global features plainly visible in the figure and typical of systems governed by nonlinear differential equations. Fortunately, laser phase diagrams can be also studied with discrete maps, not differential equations. Two examples are the Ikeda map, governing the dynamics of ring cavities containing a nonlinear dielectric medium [10], and the Hénon map, known to qualitatively reproduce the organization observed in CO 2 lasers [3]. By iterating equations recursively, discrete maps provide direct access to the whole spectrum of periodic points underlying the stability skeleton of the system. The quest for understanding the origin and organization of laser oscillations leads quickly to the necessity of understanding the exact nature of the numbers underlying the solutions of the equations of motion. As discussed below, it is the number-field infrastructure that provides the scaffolding for the physical phasespace and phase diagrams.
The exceptional periodic points defining 'centers' of stability in phase diagrams grow like weeds, seeming to obey no other law than that of chance, and nobody can predict yet where stability regions corresponding to specific laser modes will sprout. However, much can be learned by studying the selforganization of orbital points analytically in map models, i.e. exactly, with numerical approximations. This is the approach Figure 1. Stability phase diagrams illustrating the systematic accumulation of self-similar phases of periodic oscillations, for two laser models. Reprinted from [3], Copyright (2016), with permission from Elsevier. (a) Lyapunov diagram illustrating the self-organization of laser modes in a CO 2 laser with opto-electronic feedback, as a function of the feedback gain R and bias voltage B 0 . Black represents regions of periodic modes, colors represent non-periodic laser modes (chaos). (b) Isospike diagram obtained by counting number of local maxima per period of the oscillations. Non-periodic oscillations are represented in black. (c) Lyapunov diagram as in (a), but for an optically injected semiconductor laser in an optically injected semiconductor laser, where K is the intensity of the injected field and ω is the detuning frequency. taken here. We explore one classical aspect that seems to share some similarity with what is known from quantum theory. For the wide class of systems governed by polynomial equations of motion, we argue that the quantum concept of entanglement has a purely classical analogon. Entanglement is a key quantum concept studied extensively nowadays in the quest for quantum computers, machines expected to turn current binary computing devices into sluggish oddities of the past [11,12]. There is also significant work to harness the potential of entanglement in developing systems for quantum and post-quantum cryptography, and more [13][14][15].
Specifically, this paper argues that a property analogous to quantum entanglement-arithmetic entanglement-operates at the purely classical level. It arises from ambiguities inherent to the multivaluedness normally present in the solutions of classical equations of motion. In a nutshell, the classical analogon is a manifestation of unavoidable arithmetical interconnections in phase-space (see below). Individually, solutions of classical equations of motion, classical states, contain arithmetical imprints that depend on the specific branches of the multivalued functions underlying them. For two representative classical systems, we show explicitly and exactly how this works. For each system, we derive a compact representation, a carrier, embodying simultaneously complete sets of classical states of the system. Such carriers are parameterized by the multivaluedness of the individual classical states. From such carriers one may extract any classical state whenever needed by simply selecting the suitable branch corresponding to the state desired.
The first example considered is a popular one-dimensional system used as kernel of a device for performing general computations, the so-called reconfigurable chaotic computer [16]. In the second, a carrier is obtained for a considerably more sophisticated archetypal system in two-dimensions, used to describe the control space of lasers and other devices [3,17].
In both examples, selection of a single classical state from the carrier is tantamount to selecting (fixing) their duals unambiguously. The ability of rapidly and efficiently switch between classical states in nonlinear systems is essential for operating the aforementioned reconfigurable computers.
Before proceeding, we mention an interesting classical analogy of entanglement discussed by Spreeuw [18], using classical light beams, and an attempt by Briggs and Eisfeld [19] to simulate properties of entanglement with classical oscillators. Our approach, however, is rather distinct, combining current results from dynamical systems and arithmetical properties of the solutions of classical equations of motion. The theoretical framework needed here to establish orbital carriers amounts to solving cubic polynomials using elementary symmetrical functions. This modest framework works for any arbitrary dynamical systems of algebraic origin [20].
It is a great pleasure to dedicate this paper to Wolfgang Schleich. It combines topics of many of our lively and 'broad band' discussions while doing thesis work in the 'theoretical container', the barrack placed on the back of the IPP building then lodging the Projektgruppe für Laserforschung, the precursor that in 1981 became the Max-Planck Institut für Quantenoptik and that later on moved to its present building.

Multivaluedness and ambiguity
As mentioned, arithmetic entanglement originates from the innate and inextricable ambiguity present in the solutions of classical equations of motion. To grasp this, consider the elementary problem of locating in phase-space a point X whose coordinates are given by the roots of a quadratic equation, say x qx p 2 0 2 -+ = . For definiteness, suppose further that q and p are real parameters and that q 2 −p>0. Then, abbreviating and , the location of the point X is normally given as, say, X≡X(x 1 , x 2 ). While this representation is not formally wrong, the precise location of X still involves an intrinsic ambiguity that only disappears after selecting and fixing the branch of the multivalued square-root function which appears in x 1 and x 2 . In particular, depending on the branch choice, the same representations of x 1 and x 2 may be used to land on the conjugate point X X x x , . This simple example reminds us that phase-space coordinates involving multivalued functions, what they regularly do, depend on two essential things: (i) specific branch choices, and (ii) the instant of time when the branch is fixed. These choices are independent and fix coordinates unambiguously or, equivalently, with an obvious abuse of language, fix where and when the state is projected onto the phase space. Clearly, the considerations above remain true for multivalued functions of arbitrary degrees. In the remainder of the paper we show how the intrinsic multivaluedness may be applied to realistic problems, in the form of orbital carriers.

One-dimensional example
A simple but nontrivial example of arithmetic entanglement resulting from multivaluedness in solutions of the equations of motion is provided by the iterates of the quadratic (logistic) map, namely, by the iterates the equation [21][22][23][24]: where x 1 is an arbitrary initial condition. This map is used often as a metaphor for a variety of phenomena in physics. It involves the nonlinearity of lowest order in series expansions of generic equations of motion [21]. For instance, it composes the canonical quartic map a prototype for systems like, e.g., lasers, electronic circuits, chemical, biochemical and other flows [3,21]. Applications of the quadratic map are reviewed in a book commemorating its first 200 years of use [22]. Here, we study period-three orbits of the quadratic map. They are significant because of a celebrated theorem of Li and Yorke stating that period-three implies chaos [25], and a more powerful result due to Sharkovsky [26,27]. According to Dyson, the period-three theorem 'is one of the immortal gems in the literature of mathematics.' [28].

The period-three polynomial cluster
From equation (1) it follows that period-three orbits means having three numbers, say x 1 , x 2 , x 3 , repeating indefinitely as follows: Eliminating x x , 2 3 and writing x for x 1 one finds the equation of motion specific for period-three orbital points, namely: The quadratic factor in equation (3) is not interesting because, as seen from equations (2), it implies having x 1 =x 2 =x 3 , i.e.not a legitimate period-three orbit but just fixed points (i.e. period-one orbits). Effective information about periodthree is provided by E x 3 ( ) which shows manifestly that there are two period-three orbits. The next task is to split E 3 (x) into a pair of triplets of numbers corresponding to the two separated orbits. This is done in the next section, using orbital points to build the elementary symmetric functions.

Generic decomposition of orbital clusters
The pair of orbits mixed up in E x 3 ( ) may be sorted out systematically using the standard elementary symmetric functions as follows. Denote by ξ any arbitrary root of E 3 (x). Then, according to equation (1), a period-three orbit consists of the triplet ξ, f (ξ), f (2) (ξ)=f ( f (ξ)), which we use to build the symmetric functions The elementary theory of symmetric functions teaches that one may express any pair of coefficients in terms of the remaining member of the trio. We chose to express θ 3 (ξ) and θ 2 (ξ) in terms of θ 1 (ξ), the sum of the orbital points. Since These are the coefficients needed to define the first separated orbital equation: In a similar way, calling η, f (η) and f (2) (η)=f ( f (η)) the remaining triplet of roots of E 3 (x), one sees that they obey exactly the same functional relations above, namely θ 1 (η), θ 2 (η) and θ 3 (η). Clearly, this second triplet is in general distinct from θ 1 (ξ), θ 2 (ξ) and θ 3 (ξ), since they define different orbits. However, the same formal expression, equation (8), represents both orbits simultaneously. In other words, a single expression is capable of harboring simultaneously all possible period-three orbits of the system, a remarkable achievement for such a simple expression. The orbits involve two unknowns, θ 1 (ξ) and θ 1 (η), which may be determined by the same procedure used to get to equation (8). Thus, we compute the pair representing their symmetric sum and product, coefficients of the quadratic equation used to fix their numerical values.
Denoting the six roots of E 3 (x)=0 indistinctly by x 1 , x 2 , ..., x 6 , the symmetric functions of 1 q x ( ) and θ 1 (η) are found to be In equation (9), the last equality follows from the suitable coefficient in equation (4). The quantities above define a quadratic equation, say w w a 2 0 2 -+ -= , whose roots are θ 1 (ξ) and θ 1 (η), the values needed in equation (8) and θ 1 (η)=σ. As shown by equation (9), σ is simply the sum of the orbital points. Using the constraint S 3 (σ) we may eliminate a from equation (8) and obtain simpler carrier equations involving only σ. The choices in equation (12) give To check what was accomplished arithmetically, we compute the discriminants of j 1 (x) and j 2 (x): forms which bring to mind the fact that 'any third-degree polynomial p t t  Î ( ) ( ) which is irreducible over the rationals will have a cyclic Galois group if and only if the discriminant of p(t) is a square over  [29]. Thus, the above discriminants manifest clearly the advantage of using the σ-parameterization: It automatically separates orbits which have ab ovo a cyclic Galois group or, equivalently, whose Galois group is commutative. In other words, the σ-parameterization is an algorithm that allows one to effectively split orbits. Such separation is always possible for equations arising from iterated maps which, by construction, are Abelian equations [30].

Projections of classical states in the partition generating limit
After obtaining the carriers j 1 (x) and j 2 (x) in equations (13) and (14), it is natural to ask what sort of information are encoded by them. So far, a was an arbitrary parameter. For definiteness, we now focus on the so-called 'partition generating limit' a=2, the limit found appropriate to emulate logic gates, to encode numbers, and to perform specific arithmetic operations on these numbers [16,23,24].
As seen from equation (11), in the a=2 partition generating limit orbits may be reached in two ways, either by taking σ=0 or σ=1. Denoting by Φ(x) and x F( ) the macroscopic pair of period-3 orbits, they may be reached in two distinct ways as 'projections' of the microscopic carriers j 1 (x) and j 2 (x) in equations (13) and (14): These functions may be regarded as a sort of 'dynamical conjugates' despite their different-looking functional forms and distinct discriminants. We found no Tschirnhaus transformation [30] connecting them. The above polynomials are known in the literature as Gauss' cyclotomic period equations [31,32] underlying the division of a circle into nine and seven parts, respectively [33,34]. Incidentally, Cayley [34] states that 'Many of [Cayley's] results accord with, and furnish exemplifications of general theorems contained in Jacobi's memoir Über die Kreistheilung und ihre Anwendung auf die Zahlentheorie [35]'. The above expressions demonstrate something that we have not been able to find in the literature, namely that rather than being two independent entities, the above cyclotomic period equations are in fact just distinct 'projections' originating from degenerate period-three orbital carriers, either j 1 (x) or j 2 (x). Technically, although defining non-isomorphic number fields, the polynomials corresponding to equations of motion may be nevertheless dynamically conjugated. For period four we find a set of three dynamical conjugates; for period five, there are six; for period six, there are nine; etc.
In phase-space one normally works with the macroscopic objects Φ(x) and x F( ). But microscopically the description may be done equally well using either j 1 (x, σ) or j 2 (x, σ). As seen, each σ-parameterized expression encodes or, equivalently, carries the same complete set of solutions of the equations of motion. Although orbital carriers do not seem to be linear combinations of basis functions, it is nevertheless possible to extract from them the complete set of solutions. For higher periods, conjugate families may contain much more states and more complicated σ-parameterized carriers. This is a generic property of algebraic equations of motion, not a peculiarity of the illustrative example considered. By iterating polynomials one avoids the necessity of determining individual orbital points, given by the roots of the polynomials. This means representing the dynamical evolution directly with P k (x) and S k (σ), for any arbitrary period k. Noteworthy is the fact that the temporal arbitrariness in selecting a classical state by fixing a particular value of σ resembles what happens in Wheeler's delayed choice experiment.

Orbital poly-isomorphisms
To locate precisely the nucleation of stability regions it is necessary to know about number-field infrastructure underlying the equations of motion. This may be understood considering figure 2 which presents stability phase diagrams similar to the ones in figure 1, for the quartic map This two-parameter system displays the same systematic sprouting of stability observed for the Hénon map, defined below in equation (23), and other maps [7][8][9]. The quartic map is particularly useful because of the symmetries evident from figure 2(a). In this figure, black and purple regions represent the magnitude of the multiplier m of the map [21,36], a mean-field-like quantity which samples the map derivative along all points of a given trajectory. Stability boundaries for k-periodic orbits may be obtained by sweeping parameters a and b and recording the values for which m k =1 (saddle-node bifurcations) and m k =−1 (pitchfork bifurcations). When m k =0 the k-periodic orbit is called superstable [8]. Using black for multipliers lying in the interval −1m k 0 and purple when 0m k 1, we get the diagram seen in figure 2(a). White is used to represent nonperiodic (chaotic) orbits, while yellow represents parameters for which iterations diverge. For a discussion of the early days of complex dynamics, starting in 1870 with works by Schröder [36], see the nice survey by Cremer [37]. For the period between 1906 and 1942 consult [38].
The systematic of the nucleation of periodicity islands along the path a b b 2 = is illustrated in figure 2(b), a magnification of box in figure 2(a). In it, periodic islands are represented by colors and numbers indicating the period of the first few more salient islands. The regularities in figure 2(b) are typically observed in many other models, independently if governed by maps or by differential equations [3]. Behind the nucleation of stability one finds specific sequences of numbers. For instance, the stability centers, characterized by doubly superstable orbits and indicated by white dots in figure 2(b), are all defined by real numbers which obey self-inverse twin equations [9] a b a b a ... , where the number of radicals taken in each island is the same -. White dots mark the doubly superstable points at the center of the islands [7,8]. See text. and is related in a simple manner to the orbital period k. For example, the stability center for period-one orbits is given by (a, b)=(0, 0), the intersection point of the self-reciprocal parabolas a b =  and b a =  or, equivalently, b=a 2 and a=b 2 , clearly visible in figure 2(a). Twin number satisfying equations (18) and (19) regulate the nucleation of all periodicity islands of the quartic map [8,9], not only those along the path a b b 2 = shown in the figure. On both sides of the large green period-three shrimpnearly at the center of figure 2(b) one finds a pair of periodfive shrimps. Is it possible to find transformations allowing one to move from one period-five center to the other one? Investigation of this possibility leads to a new and interesting class of problems connecting dynamics and numbers which, for simplicity, we illustrate briefly in the context of the logistic map, equation (1). The points of the simplest periodfive orbit of the map are roots of the quintic  10  42  39  5  5  28  39  10  2  7  13  5  2  1,  5  8  10  32  19  10  48  57  4  5  32  57  8  32  8  19  4  32  23,  5  13  10  52  61  10  78  183  119  5  52  183  238  70  13  61  119  70 23,  -+  -+  --+  -+  -+  ---+ - where n is an arbitrary number, positive or negative. Many additional sequences of polynomials exist, all sharing the same discriminant Δ V =11 4 . Now, back to the question of transformations among stability regions, one sees that even for the simple transformations in equation (22), the number of quintics sharing the discriminant Δ V is too big to be tested individually, specially since the number of quintics grows fast when the degree of the equation of motion grows. For the logistic map, in addition to Vandermonde's celebrated quintic in equation (20), there are five additional period-five orbits. But they turn out to be mixed up in polynomials of degrees 10 and 15, meaning that quintics corresponding to individual orbits will not have integer coefficients. For the quartic map of figure 2 the situation is much worse: there are 204 period five orbits, all mixed together in a polynomial of degree 1020, that current commercial software is not yet able to handle for arbitrary parameters a and b, or along the line a b b 2 = -. However, the fact that this polynomial cluster does not break up looks promising, indicating that it is a cluster of 204 quintics sharing the same discriminant, and eventually allowing transformations among them. We are currently developing ad hoc software to address this and related questions. They bring number theory into the realm of physics, into practical applications for determining the stability of accumulation sequences of periodic orbits.

Two dimensional example
Is it possible to obtain orbital carriers for higher-dimensional systems? Certainly: in this case, the iterative procedure results in considerably more intricate polynomial clusters but, as before, it is possible to derive exact σ-parameterized solutions. To show this, we consider the so-called Hénon map, a paradigmatic textbook example introduced originally as a model of galaxies but whose phase diagram was found subsequently to represent well the stability diagram of certain class B CO 2 lasers [3,17], and other practical devices [21]. The Hénon map is a dissipative system defined by the Parameterized equations covering the Hamiltonian limit (b 1 = -) valid for every period up to 20 are known [39]. As expected, in the b=0 limit the above equations correctly reduce to the results found previously, for the one-dimensional example. Thus, adding more parameters and/or extra dimensions only alters coefficients, not substance. Of course, higher-dimensional systems demand longer algebraic work but should present no intrinsic difficulties.

Conclusions and prospects
The above examples provide explicit evidence that an analogon of the quantum concept of entanglement may be identified in purely classical states, thanks to the intrinsic multivaluedness prevalent in the solutions of the equations of motion. Multivaluedness was used to extract convenient orbital carriers subsuming the set of possible states into a single expression from which every state possible may be conveniently recovered when desired. The carriers in equations (13) and (14) revealed a startling result: pairs of the celebrated Gauss' cyclotomic period equations [31,32] are nothing but just distinct projections derived from more fundamental entities, j 1 (x) or j 2 (x), as indicated plainly in equations (15) and (16). Transformation properties of orbits sharing the same discriminant were also discussed. Although promising research avenues exist, they are not yet easily accessible with standard software for algebraic manipulation currently available but might be within reach with specially designed programming.
In some sense, carriers like x 1 j ( ) and j 2 (x) play a role akin to quantum linear combinations and their associated nonlocality features. These carriers offer a fresh view concerning the origin of Einstein's spukhafte Fernwirkung, the 'spooky action at a distance', an effect demonstrated repeatedly in experiments but not yet fully understood theoretically [11]. A recent paper reported compelling evidence that entanglement is an inescapable feature of any physical theory aspiring to be more general than classical physics, and that viable candidates of extended theories must contain classical physics in a suitable limit [40]. In contrast, the approach envisioned here seems to indicate that, rather than arising from a top-down design, entanglement might have been always lurking at the classical level, disguised in the form of branch ambiguities, as demonstrated here explicitly for algebraic systems.
The proof-of-principle that carriers can mimic some aspects of quantum physics opens a few avenues for fruitful applications in classical physics based on considerations of full orbits, not sequences of orbital points as it is normally done. The analogy of carriers and quantum phenomena discussed here is admittedly no more than an analogy, whether a deep or a specious one remaining to be further explored. However, the important point is that, totally independent of analogies, orbital carriers have a solid foundation on the theory of equations and are ready to be used in other systems. Carriers provide an efficient technique for tracking and controlling orbits of generic systems governed by polynomial equations of motion. We believe that many surprising results based on them are awaiting discovery.