Chaotic dynamics in spinor Bose–Einstein condensates

Concentrating on experimentally relevant cases we discuss nonlinear dynamics and deterministic chaos in F=1 and F=2 spinor Bose–Einstein condensates. A thorough numerical and analytical treatment of the equations of motion of symmetric spin dynamics in the single-mode approximation elucidates the complex dynamics and uncovers regions of classical chaos in this important quantum gas system.


2
Since the groundbreaking demonstrations of Bose-Einstein condensates with spin degree of freedom as quasi spin- 1 2 [1], spin-1 [2] and spin-2 systems [3] and their theoretical analysis [4]- [8], spinor condensates have become a focal point of modern ultracold quantum gas research [9]- [21]. The possibility of tailoring all internal and external parameters has made these systems near ideal experimental quantum magnetic model systems, covering the full range from magnetic phases to spin correlations. Fascinating examples are the realization of quantum control of spinor dynamics [22] and the discovery of spontaneous fragmentation [13] as well as a magnetically tuned resonance [18]. Controlling the temperature of the system, even multicomponent thermodynamic [16,23] and thermal decoherence effects [17] become accessible in experiments. These fascinating experimental possibilities are complemented by a detailed theoretical understanding of F = 1 2 and 1 dynamics, which can be based on relatively simple principles [24]. In many cases the dynamics can even be accessed via analytical expressions [11,17]; however, the dynamics at higher spin values becomes increasingly complex and less intuitive.
In this paper, we show that going from F = 1 to F = 2 systems indeed adds nontrivial and perhaps surprising new aspects to spinor dynamics, such as the occurrence of classical chaos, which might have significant consequences for the interpretation of experiments and the development of spinor correlations. We discuss the phenomenon of chaotic spin dynamics using the initially fully stretched state oriented perpendicular to the magnetic field as an example analogous to the experiments in [17,18]. Working within the single mode approximation we identify different parameter regimes based on numerical simulations of the 'classical' meanfield spinor equations.
Spinor condensates add new aspects to the large field of 'spin dynamics', with relations mainly in the context of classical and quantum chaos [25,26], nuclear magnetic resonance (NMR) [27] and solid-state magnetism [28]. On the other hand, spinor Bose-Einstein condensates are closely related to single-component Bose-Einstein condensates in multiplewell potentials [29].
The paper is organized as follows: in section 1, we show that neglecting higher order effects one can introduce general spinor mean field equations for an arbitrary spin. Within this approximation we discuss the similarities and differences between the F = 1 and F = 2 cases as particular examples. In section 2, we introduce the magnetic-field-dependent single-mode spin dynamics in these two cases and show that the F = 2 system exhibits regions of surprising irregularity, which as a matter of fact also show a nonvanishing Lyapunov exponent. Section 4 is dedicated to the more detailed discussion of these regions, using appropriate Poincaré sections of the phase space, which clearly demonstrate the interruption of chaotic dynamics by 'islands' of regularity.

Mean field theory of spinor condensates
In order to introduce the notation and formally derive the mean-field equations for spinor condensates, it is instructive to first review the analogue in the single-component Bose-Einstein condensate theory. The Bose-Einstein condensates in dilute gases of weakly interacting particles are typically described by the Gross-Pitaevskii equation (see e.g. the review article [30] and references therein), 3 In this version for a single-component or spin-0 gas, ( r , t) is the space-and time-dependent order parameter of the condensate, which can be interpreted as a wavefunction or state vector in the single-particle Hilbert space. The Gross-Pitaevskii equation (1) is then identified as the single-particle Schrödinger equation with an additional nonlinear term representing the average effect of all remaining particles on any single one (mean field). Formally, the Gross-Pitaevskii equation is derived from the energy functional as the variational derivative with respect to * . Note that is normalized to the total number of particles N , thus | | 2 is the density of particles. The interaction of particles, described by a two-body potential V int ( r 1 − r 2 ), is reduced to an effective contact interactionhgδ( r 1 − r 2 ) in the s-wave scattering limit, applicable in a dilute gas at low energies. The coupling constant g = 4πha/m is related 4 to the s-wave scattering length a.
Generalizing to spin-F particles, ( r , t) is replaced by a (2F + 1)-component wavefunction m ( r , t), m = −F · · · + F. If an additional external potential depends on the spin orientation, it is also replaced with a 2F + 1-component version.
Owing to rotational symmetry and the necessity for the pair wavefunction of the colliding bosons to have even parity, the scattering length depends only on the total spin F tot = 0, 2, . . . , 2F of a pair of colliding particles [4], i.e. there are F + 1 potentially different scattering lengths a 0 , . . . , a 2F , and the effective interaction potential takes the form where P f are projectors on the subspace F tot = f . It is useful to rewrite the projectors as linear combinations of the identity operator, F 1 · F 2 and P 0 [4,5,7], where the coupling constants g n depend on the scattering lengths a f (g 2 = 0 for the case of F = 1) as given in appendix B. Using this interaction potential for a Bose-Einstein condensate with spin degree of freedom, the energy functional becomes where The multicomponent Gross-Pitaevskii equation for spin F = 1 and F = 2 is easily (though lengthily) recovered by taking the derivative of this energy functional ih˙ k ( r , t) = ∂ H /∂ * k ( r , t), leading to a system of 2F + 1 coupled nonlinear partial differential equations. However, the resulting general equations are rather complicated and do not yet include experimentally important magnetic field effects. In order to gain an intuitive insight into spinor dynamics in magnetic fields, we will in the following introduce some approximations and specializations before considering Zeeman contributions to the energy due to an external magnetic field.
As a first step, we concentrate on the intrinsic spinor dynamics and neglect spatial variations, i.e. we adopt the assumption that spatial and spin degrees of freedom decouple. This is expressed in a product ansatz for the wavefunction, This ansatz is known as the single-mode approximation [6] (also used implicitly in [4,5]), and is characterized by the assumptions that φ is fixed in time and ζ k does not vary in space. This approximation is generally considered appropriate if the spin-dependent part of the interaction energy and the particle number are small, i.e. g 1,2 g 0 and N < 10 4 [6,31] (also compare [32] for a more detailed analysis of the validity of the SMA). Summing over m F states and integrating over space, respectively, the energy functional (5) decomposes into a spin-independent part that takes the same form as for a single-component gas, and a spin-dependent part that takes the average density n ≡ |φ| 4 d 3 r/N as a parameter, The effect of a magnetic field B 0 in the z-direction can now be added in the form of the linear and quadratic Zeeman effects, which enter into (11) as an additional external potential which is homogeneous in space but depends on spin orientation, The parameters p and q, with | p| = (1/2)(µ B /h)B 0 and |q| = p 2 /ω hfs , are both positive (negative) for F = 1 (F = 2). The hyperfine splitting is ω hfs ≈ 2π × 6.824 GHz for 87 Rb. The equations of motion resulting from this energy functional conserve the norm m |ζ m | 2 , the energy H spin and the z-component of the spin F z = m m|ζ m | 2 . As a result of spin conservation, the linear Zeeman effect is irrelevant for dynamics and can be neglected in the spin-dependent energy functional. The equations are then also invariant under the 5 transformation T : ζ +m → ζ −m , observing that for m → −m, F z,y → −F z,y and F x → F x as well as S 0 → S 0 . T corresponds to a half-cycle rotation of the coordinate system around the x-axis, T ∝ exp(iπF x ). Obviously, T = T −1 is its own inverse and has eigenvalues ±1.
In order to take advantage of this symmetry, we will concentrate on the dynamics within the symmetric subspace defined by T|ζ = |ζ , implying that F z = F y = 0. As a last approximation we further neglect g 2 terms in F = 2 dynamics, which allows us to cast the spindependent energy functional in the general form 5 This equation is the basis for the discussion and simulations in the following sections. It is remarkable that despite its relatively simple form, it already contains the essential step to chaotic behavior, when going from F = 1 to F = 2 systems. Since (14) as well as the equations of motion for arbitrary observables A (which can be derived from (14) noting that˙ A = ζ |A|ζ + ζ |A|ζ ), are formally independent of the value of F, the question arises how differences between F = 1 and F = 2 come into play. The answer lies in the fact that the equation of motion for any observable will generally contain other observables which are products of the operators involved, and will thus not lead to a closed system of equations. As an example that will also be of use later, we consider F x : The equation of motion of F z F y + F y F z can be calculated along the same lines, but will contain products of three operators, and so on. However, since we are operating in a finite-dimensional Hilbert space, all observables can be expressed as linear combinations of a finite number of basis operators in this Hilbert space. Doing this, one arrives at a closed set of at most (2F + 1) 2 equations of motion. Symmetry with respect to m → −m further reduces the number of basis matrices to (F + 1)(2F + 1). Conservation of energy and particle number (or normalization), the freedom of choice of a global phase, as well as the fact that only spin operators and their products appear in the equations of motion, lead to differential equations with only two (in the case of F = 1) and four (F = 2) degrees of freedom. Thus, chaotic dynamics becomes possible (and is indeed realized as we will discuss in the following) for F 2 but not for F = 1 [33].

Single-mode dynamics of spinor condensates
We now turn to the discussion of the dynamics resulting from (14) in the cases F = 1 and F = 2, highlighting common characteristics as well as distinctive features arising from the larger number of degrees of freedom in F = 2.
6 Dynamics in single-mode approximation is equivalent to the homogeneous spinor Bose gas if we assume that φ( r ) remains fixed, i.e. in the ground state of H 0 . 6 The time-dependent Gross-Pitaevskii equation for the spin components ζ k is obtained from (14) by differentiation with respect to ζ * k , Using H = H spin /h N , the equation of motion for the spinor |ζ becomes or more explicitly for F = 1 with the abbreviations . As discussed in detail in [17], this dynamics is characterized by the competition between interaction energy, parameterized by g 1 , and the (quadratic) Zeeman energy q F 2 z (see also (14)). The parameter regimes where one of the two dominates will be identified in the following as the interaction regime (|g 1 n | |q|) and the Zeeman regime (|g 1 n | |q|), respectively.
The solution of the equations of motion (19) and (20) is in general rather complicated due to their nonlinearity. In order to gain insights into the behavior of the system, we will therefore concentrate on a special case, for which analytical results can be obtained. For a number of reasons, the fully transversely magnetized state |ζ π/2 = exp(iπF y /2)| − F , belonging to the class of symmetric states, is of particular interest both experimentally and theoretically. Firstly, it is easily and reliably prepared and has been used in a number of experiments [17,18]. Secondly, it is a state of maximum interaction energy for anti-ferromagnetic species (e.g. sodium in F = 1 [34] and rubidium in F = 2 [35]) and as such is expected to exhibit rich dynamics. Note that, as we will see in section 4, small experimental deviations from this special initial state  will not modify the general chaotic behavior of the system and thus will not impede possible observation of these effects. For F = 1, conservation of spin and energy as well as normalization and global phase invariance reduce the number of degrees of freedom to just two, e.g. the population ρ o of m F = 0 and a relative phase θ ≡ θ +1 + θ −1 − 2θ 0 [11], where ζ m = √ ρ m e iθ m . This system is integrable, and in fact for the particular initial condition ζ (0) = ζ π/2 = (1/4, 1/2, 1/4) T , an analytic solution can be obtained [17,36]. The time-dependent populations of the spin components can be written in terms of Jacobi elliptic functions [37] |ζ where k = g 1 n /q is the ratio of interaction energy to the quadratic Zeeman effect. Figure 1 illustrates the analytic solution across a wide range of k and in the particular case k = 1. For small |k| 1 (the Zeeman regime), the Jacobi elliptic functions can be approximated by ordinary trigonometric ones, i.e. sn k (x) ≈ sin(x), cn k (x) ≈ cos(x), dn k (x) ≈ 1. Thus, the populations of (22) oscillate with an amplitude given by k, and a period which is π/q for small |k|. In the opposite regime, |k| 1 (interaction regime), the identity sn k (x) = 1 k sn 1/k (kx) leads to another trigonometric approximation, sn k (x) ≈ sin(kx)/k. In this regime, populations oscillate with amplitude 1/k and period π/(g 1 n ). The crossover region |k| ≈ 1 (crossover regime) exhibits a maximum of the amplitude, while the oscillation period diverges. In other words, at |k| = 1 the evolution becomes aperiodic, and the m F = 0 population asymptotically approaches unity.
Owing to the larger number of degrees of freedom, the equations of motion for 87 Rb F = 2 are not integrable and no exact analytic solution exists. However, in the limiting cases of very small and very large quadratic Zeeman effects q, approximate solutions for the initial condition ζ (0) = ζ π/2 = (1/4, 1/2, √ 3/8, 1/2, 1/4) T have been obtained [18,36] as well, exhibiting similar oscillatory behavior as in F = 1. In the Zeeman regime (|q| → ∞), the fundamental frequency of oscillation is q/π with an amplitude ∝ g 1 n /q. In the opposite limit, the interaction regime |q| → 0, the frequency of oscillation 2g 1 n /π is determined by the same mean-field interaction term but is twice as fast compared to F = 1.
By analogy to the F = 1 case, and considering the limiting cases, a resonance-like phenomenon in the crossover region is also expected in F = 2 and indeed has been found in Step-wise minima Step-wise maxima Single run 0 100 200 300 400 500 600 0 Figure 2. Simulated F = 2 population dynamics at g 1 n = 50 s −1 and q = 2.5g 1 n ; the initial state ζ π/2 is the fully transversely magnetized state. Note how the indicated 'step-wise' maxima and minima and the 'single-run' trajectory start to diverge as a result of accumulated numerical errors. The 'step-wise' and the 'single-run' trajectory are calculated using the same code with a different choice of interpolation points.
experiments [18]. However, a numerical simulation of a typical trajectory (figure 2) from this region, starting from the fully transversely magnetized state, is far from the regular behavior in the F = 1 case. It instead shows a strongly irregular oscillation with ill-defined amplitude and frequency. The calculated trajectory is also extraordinarily sensitive to unavoidable numerical errors, which is a typical sign of deterministic chaos. This is a first indication of nontrivial differences between F = 1 and F = 2 dynamics, which we will further elucidate in the following. The exponential amplification of small deviations is a fundamental characteristic of chaotic systems and necessarily introduces an element of randomness into numerical studies of such systems. Therefore, statistical tools such as histograms and averaged quantities (often assuming ergodicity) are routinely employed in their analysis.

Chaotic spin dynamics
In order to further analyze the magnetic field-dependent nature of the spinor oscillations in terms of 'irregularity', we have performed numerical simulations extracting amplitude and cycle time. The analysis is performed in analogy to obtaining the amplitude and cycle time from welldefined periodic signals by taking the differences in value and time of each pair of minimum and maximum. Applying this method of measurement to an irregular signal such as the m F = 0 population trajectory ρ 0 as in figure 2, we obtain a different amplitude and period for each minimum-maximum pair 7 : where ρ 0 (t) has a maximum at t = t max , equivalent to its time derivative having a zero crossing from positive to negative,ρ 0 (t max ) = 0 (t min correspondingly).
Analyzing trajectories over a given number of oscillations (more than 3000 minimummaximum pairs) we obtain a histogram of amplitudes and periods as in figure 3. The plot shows, for a range of given parameters |q| (the horizontal axis), the frequency of occurrence of an amplitude or period of a given value (the vertical axis) as color-coded pixels. For these simulations, |g 1 n | = 1 has been assumed, which merely fixes the timescale. First of all, figure 3 demonstrates once again the prevalence of simple oscillations of decreasing amplitude in the limiting cases |q| → 0 (the interaction regime) and |q| → ∞ (the Zeeman regime). In order to classify regions in the histograms as chaotic or regular, the graphs also show the Lyapunov exponent of the trajectories, which has been calculated according to appendix A. Chaotic trajectories, as opposed to regular ones, are characterized by a nonzero positive Lyapunov exponent; numerically however, the line has to be drawn at a small finite value. In figure 3, light red shaded parameter regions have a numerically calculated Lyapunov exponent of 0.01 or larger.
The additional degrees of freedom in F = 2 cause a broadening of the histogram due to pseudoperiodic motion, in particular in the wings towards the Zeeman and interaction regimes, or a completely smeared-out distribution in the chaotic regimes. This smearing-out is a common characteristic of the chaotic regime, as well as the interspersed islands of order, where the motion becomes (pseudo-)periodic for a narrow range of parameters q. One of them, at |q| = 2.82, will be studied in more detail in the following using a Poincaré surface of section to visualize trajectories.
Fitting power laws A(q) = aq α and T (q) = tq β to the center-of-mass of the wings of the histograms, we recover the asymptotic behavior known from the analytic solution in F = 1, and expected from perturbative solutions in F = 2 [18]. The oscillation period is constant and of the order of one (or g −1 in general) in the interaction regime, and follows q −1 in the Zeeman regime. The oscillation amplitude is maximal in the crossover region and decays both in the Zeeman regime and in the interaction regime. In the latter case, the best-fit exponent α ≈ 1 is in quantitative agreement with the perturbative solutions. In the Zeeman regime, the best-fit exponent deviates from the expected value α = −1, probably because at |q| = 10 the limiting case is not yet fully realized.
In summary, the histograms of F = 2 spin dynamics starting from the fully transversely magnetized state reveal both regular oscillatory behavior in the wings and a broad distribution of amplitudes and periodic times in the crossover region. While the former is fully analogous to F = 1, the latter is a characteristic feature of F = 2 and supports the existence of deterministc chaos in this system. The following section reveals that irregularity and positive Lyapunov exponents are not limited to the particular initial condition of the fully transversely magnetized state, but are indeed a feature of the whole system in certain parameter ranges.

Phase space and Poincaré sections
In order to gain further insights into the chaotic parameter region, the 'robustness' of chaos and in particular the above sequence of chaotic and regular ranges, it is useful to define and analyze suitable Poincaré sections. Poincaré sections have proved a powerful tool for the visualization of   complex dynamics in phase space. A Poincaré section reduces a time-continuous n-dimensional system to a discrete n − 1-dimensional one by looking at the intersections of trajectories with a given surface in phase space [38]. The resulting discrete dynamical system still contains the essential features of the continuous one. The main task is to pick suitable observables that allow us to define a Poincaré section of physical relevance. In the case of a 3D phase space as in the present case, the surface of section is 2D and can be easily visualized.
From the perspective of dynamic systems, using observables instead of spinor components is just a change of variables, and it is possible to fully describe the system (20) in terms of six suitably chosen observables instead of the real and imaginary parts of three spinor components ζ n = ζ n + iζ n , n = 0, 1 and 2. Naturally, one chooses known constants of the motion as new variables. These are the 'number of particles' N and the energy E, Additionally, we choose F x ≡ F x and F zy ≡ F z F y + F y F z /2 due to their intuitive meaning as the transverse spin and its time derivative, Finally, for compatibility with previous results, we take the population and phase of the m F = 0 component, θ 0 is the only quantity influenced by the choice of the global phase, and consequently is arbitrary. Ignoring θ 0 amounts to projecting the phase space onto a plane of constant θ 0 , and the freedom of choice of the global phase ensures that the physically meaningful dynamics is fully captured in the projection. Having found suitable variables and a Poincaré surface of section, we can apply this technique to gain a deeper knowledge of the dynamical structure of our system. First of all we note that the Poincaré surface of section F zy = 0 corresponds to points of the trajectory where F x is minimal or maximal (turning points).
for any trajectory (F x (t), ρ(t)) there is a mirror image (−F x (t), ρ(t)). Minima and maxima of F x are reversed in the mirror image. Thus, instead of plotting a Poincaré section of e.g. minima (zero crossings of F zy with positive slope) in the full range F x = −2 . . . 2, it is equivalent to plot only half the range but using both minima and maxima. Also, F 2 z = (g 1 n F x 2 /2 − E)/|q| may be used instead of F x when the sign is not important. In the following, we apply both ways of plotting the Poincaré section to numerical solutions of (20) in a series of parameter values |q|, while g 1 n = 1 without loss of generality 8 . We look at the energy shell E = 2g 1 n − |q|, which is the energy of the fully transversely magnetized state. In this case, the relation of F 2 z and F x is We obtain trajectories by numerically solving (20), using a variable time step method and intercepting the zero crossings of 1 2 F y F z + F z F y . For a given initial state, we thus obtain a series of 1000 intersections of the trajectory with the corresponding Poincaré surface of section. The spinor amplitudes at the intersection points are then transformed to (ρ 0 , F x ) phase space. Initial conditions for the trajectories are chosen on a regular grid in (ρ 0 , F x ) space, limited to the valid regions indicated in figure C.1. Additionally, a number of randomly chosen trajectories starting from the neighborhood of the fully transversely magnetized state are traced (referred to as 'neighbourhood trajectories' in the following). The initial distances of the starting points from the fully transversely magnetized state cover four orders of magnitude, thus rendering a complete picture of the phase-space structures the fully transversely magnetized state might belong to.
For visualization, each intersection of a trajectory with the Poincaré surface of section is plotted as a small gray dot in the (ρ, F qz ) plane, its shade corresponding to the Lyapunov exponent of the whole trajectory it belongs to (figure 4). Intersections of neighborhood trajectories are plotted as larger dots colored red or green according to the direction of sign change of the zero crossing of 1 2 F y F z + F z F y or, equivalently, depending on whether the intersection constitutes a maximum or a minimum of |F x |. Figure 4 presents three examples of such phase-space portraits, elucidating the technique and its physical interpetration. At |q| = 1.4, well below the crossover region in the interaction regime (compare figure 3), phase space 9 is filled with pseudoperdiodic trajectories, showing up as circular structures on the surface of section. The Lyapunov exponent is zero throughout the whole phase space, indicating stable trajectories. The fully transversely magnetized state is part of a pseudoperiodic cycle, and all neighborhood trajectories trace out the same structure, confirming its stability. Thus, small deviations in the preparation of the fully transversely 13 magnetized state will not affect its characteristic dynamics. In contrast, in the crossover region at |q| = 2.3 and |q| = 2.9, large parts of the phase space are densely filled with trajectories having positive Lyapunov exponents. At the same time, the strong mixing associated with a positive Lyapunov exponent washes out any visible structure in these parts of the phase space. Other parts remain regular-at the same parameter q, pseudoperiodic and chaotic trajectories coexist, depending on the initial condition.
The two phase-space portraits from the crossover regime also demonstrate a characteristic property of conservative chaotic systems: regular and chaotic regions are densely intermingled at all scales. As an example, at |q| = 2.3, the fully transversely magnetized state is part of a chaotic region, and its instability shows up in the spreading of neighborhood trajectories throughout the chaotic part of phase space. In contrast, at q = 2.9, the fully transversely magnetized state happens to be part of a small 'island of order' hidden between mostly chaotic trajectories. Figure 5 is a close-up of this particular 'island of order' across a parameter range of |q| = 2.82 . . . 2.86. Looking at the whole phase space (figure 5, middle column), we notice that its overall structure does not change over the narrow range of parameters chosen here. There is a large part densely filled with chaotic trajectories, changing shape only very slightly. However, the behavior of neighborhood trajectories varies dramatically with q, switching from chaotic motion to near-periodic cycles. This change in behavior shows up in the histograms (figure 3) as a sudden collapse of the broad amplitude and period distributions to only a few discrete values. The reason only becomes obvious in a close-up of the neighborhood of the fully transversely magnetized state (figure 5, right column): trajectories become trapped when the fully transversely magnetized state happens to be close to a periodic cycle, which passes by as q varies. At |q| = 2.850, the fully transversely magnetized state nearly coincides with the intersection of the periodic cycle and the surface of section.
The fully transversely magnetized state is an intersection point of maximal |F x | = 2, but exactly the same behavior is seen for the corresponding turning point of minimal |F x | (figure 5, left column). Interestingly, part of the structure surrounding the periodic cycle seems to extend across the border of the physical phase-space region. This is seen particularly clearly in the circular path traced by intersections of maximum |F x | at |q| = 2.86, half of which is at the high-F x boundary and half of it at the low-F x boundary of the permitted region. This points to a more complicated topology of the accessible phase space than is obvious in this visualization. There are points at the low-F x and high-F x boundaries that must be dynamically equivalent.

Conclusion and outlook
In conclusion, we have demonstrated the emergence of classical chaos in F = 2 spinor dynamics. In particular, we have classified the corresponding parameter ranges for an experimentally relevant case (the fully transversely magnetized state) and have analyzed the occurrence of 'islands of order' as a function of the mean field to the quadratic Zeeman energy ratio. Expanding our perspective to all possible initial conditions, we have shown that deterministic chaos is indeed a characteristic of the F = 2 system as a whole and thus robust with respect to experimental inaccuracies.
Concerning experimental observation of the chaotic regime, several issues have to be taken into account. Firstly, it is important to note that the characterization of a trajectory as chaotic or regular relies on a statistical analysis. In principle, a large number of oscillatory cycles for  a given set of parameters need to be analyzed in order to determine Lyapunov exponents or to obtain phase-space images, which certainly imposes tremendous requirements on stability and reproducibility. However, figure 2 demonstrates that already 10-20 cycles suffice to observe clear irregularities. Secondly, the influence of possible multimode dynamics in experiments has to be carefully considered. A restriction to single-mode dynamics would require spinor condensates in a trapping volume much smaller than the spin healing length ξ =h/ √ 2mg 1 n , which could be realized with low atom number dipole traps or in a large-scale optical lattice. A first step towards an experimental realization of chaotic spin dynamics might thus be a straightforward combination of existing experimental techniques: small-scale optical traps and a spin-sensitive nondestructive detection technique [12], which allows the observation of a series of oscillations for one single condensate realization (thus removing parameter uncertainties in repeated measurements with destructive imaging techniques).
Our findings provide important insights into the nontrivial dynamics of high spin spinor condensates and help to identify regimes of 'deterministic' versus 'chaotic' nature. An interesting future question is how the predicted classical chaos is modified by the quantum nature of the spinor system, which possibly opens a new avenue to experimental investigations of quantum chaos. A straightforward extension to spatially extended 'multimode' systems would then create fascinating perspectives for the analysis of the interplay of correlations and chaos in tailored spinor condensate systems.
We obtain an equation of motion for the perturbation by linearizing F in the neighborhood of the trajectory z(t), where ∂F is the Jacobian matrix of F and depends on z. The equation of motion for the perturbation is linear, but time-dependent via the trajectory z(t), Solving this equation of motion, we obtain the time evolution operator J(z 0 , t) that maps any initial perturbation ζ 0 to its value at time t, The Lyapunov exponent is a measure of how fast two infinitesimally close trajectories separate. Identifying one trajectory as z(t) and one as z(t) + ζ (t), we define the Lyapunov exponent λ λ = lim Its value can be calculated from the linear approximation (following [39]), where, in the latter variant, the length of the initial perturbation has dropped out and the expression depends only on its orientation n = ζ 0 /| ζ 0 |. The matrix M = J T J is symmetric and positive, and therefore diagonalizable with real positive eigenvalues 2 1 , . . . , 2 N . Choosing the respective eigenvectors as our initial orientation n, we obtain a spectrum of Lyapunov exponents, Whether a trajectory is stable or not is determined by the largest Lyapunov exponent λ 0 = max{λ i }. If the orientation of an initial perturbation is not exactly perpendicular to the eigenvector u 0 corresponding to λ 0 , i.e. if ζ 0 has a component in the direction u 0 , this component will grow exponentially at a rate larger than any other component of ζ 0 , and will thus dominate in (A.7). In other words, an arbitrary initial orientation (not perpendicular to u 0 ) will be tilted over time towards the direction of the maximal growth rate.
The largest Lyapunov exponent can be numerically determined by solving the equations of motion (A.1) and (A.4) simultaneously for a set of perturbations of length | ζ | = 1 and different initial orientations. Since the length | ζ | of the perturbation grows exponentially in time, it is necessary to rescale ζ to unity at regular intervals and record its logarithm. The sum of the logarithms, divided by the total time of evolution, is an approximation for the largest Lyapunov exponent. If the algorithm works as expected and produces a good approximation, the resulting Lyapunov exponent should be the same for all initial orientations (except for one, possibly).

Appendix B. Review of the mean-field description of spinor Bose-Einstein condensates
In order to calculate the spin-dependent mean-field potential, it is useful to rewrite the projectors as products of single-particle spin operators. For F = 1 and 2, the commonly used decompositions are [4,5,7] F = 1 : In the case of F = 2, the projection on total spin F tot = 0 is needed besides the identity operator and F 1 · F 2 . The interaction potential then takes the form 4πh m 7a 0 − 10a 2 + 3a 4 7 . (B.5)
which can be numerically solved for the phase differences θ 0 − θ 1 and θ 1 − θ 2 . The global phase is irrelevant, but defined by tan θ 0 for formal completeness. From the populations and phases, the complex spinor components ζ m = √ ρ m e iθ m can be recovered. The relations between this choice of quantities and the spinor components are in fact invertible, proving the change of variables. However, only a limited domain of values (N , E, F x , F zy , ρ 0 and θ 0 ) is physically accessible, and only within this domain the new variables can be traced back to spinor components. In contrast, any spinor |ζ can be converted into the new variables (N , E, F x , F zy , ρ 0 and θ 0 ).
Of the six variables (N , E, F x , F zy , ρ 0 and θ 0 ), only F x , F zy and ρ 0 are dynamically relevant. In the following, we will further eliminate F zy by defining a Poincaré surface of section F zy = 0. Still, not the full domain (ρ 0 , F x ) ∈ [0 . . . 1] × [ − 2 . . . 2] is physically accessible. Apart from the implicit relationship of F x = F x and ρ 0 (C.4), the most important restrictions result from 0 F 2 z 4 and 0 ρ 1,2 1. For the specific energy shell of the fully transversely magnetized state., E = 2g 1 n − |q| with N = 1, and g 1 n , |q| > 0, we have the following nontrivial constraints: The two conditions on ρ 0 coincide only in the limiting case |F x | = √ 4 − 2|q|/g 1 n , where ρ 0 = 1 necessarily. For |q|/(g 1 n ) > 2, the conditions remain valid but do no longer constrain F x . Figure C.1 illustrates the domains of validity for the different cases.
On the Poincaré surface of section F zy = 0, even further restrictions apply as indicated in figure C.1. The shaded areas in figure C.1 have been calculated numerically by trying to invert (29) on a dense grid in (ρ 0 , F x ) phase space.
Note the different meanings of allowed domains on the Poincaré section, as opposed to a projection of the full phase space. In the former case, a trajectory may 'jump' from one domain to another, even if they are not connected. However, the projection of the trajectory between intersections with the Poincaré surface is still subject to conditions as explained previously. Thus for |q|/g 1 n < 2, F x can never change sign, but for |q|/g 1 n > 2 it does even though the dark shaded areas in figure C.1 may not actually meet.