Boundary time crystals in collective d-level systems
Luis Fernando dos Prazeres, Leonardo da Silva Souza, Fernando Iemini
BBoundary time crystals in collective d -level systems Luis Fernando dos Prazeres, Leonardo da Silva Souza, and Fernando Iemini Instituto de F´ısica, Universidade Federal Fluminense, Niter´oi, Rio de Janeiro, Brazil (Dated: February 9, 2021)Boundary time crystals (BTC’s) are non-equilibrium phases of matter occurring in quantumsystems in contact to an environment, for which a macroscopic fraction of the many body systembreaks time translation symmetry. We study BTC’s in collective d -level systems, focusing in thecases with d = 2, 3 and 4. We find that BTC’s appear in different forms for the different cases. Wefirst consider the model with collective d = 2-level systems [presented in Ref.[1]], whose dynamics isdescribed by a Lindblad master equation, and perform a throughout analysis of its phase diagramand Jacobian stability for different interacting terms in the coherent Hamiltonian. In particular,using perturbation theory for general (non Hermitian) matrices we obtain analytically how a specific Z symmetry breaking Hamiltonian term destroys the BTC phase in the model. Based on theseresults we define a d = 4 model composed of a pair of collective 2-level systems interacting witheach other. We show that this model support richer dynamical phases, ranging from limit-cycles,period-doubling bifurcations and a route to chaotic dynamics. The BTC phase is more robust inthis case, not annihilated by the former symmetry breaking Hamiltonian terms. The model withcollective d = 3-level systems is defined similarly, as competing pairs of levels, but sharing a commoncollective level. The dynamics can deviate significantly from the previous cases, supporting phaseswith the coexistence of multiple limit-cycles, closed orbits and a full degeneracy of zero Lyapunovexponents. I. INTRODUCTION
The classification of different phases of matter accord-ing to their spontaneous symmetry breaking (SSB) is acornerstone of physics and one of Landau’s legacy [2, 3].It is based on the idea that the system, in the thermo-dynamic limit, can break some of its symmetries due tothermal or quantum fluctuations giving rise to differentphases of matter, as e.g. crystals in case a spatial transla-tional symmetry is broken, superfluids for gauge symme-tries, ferromagnets in the case of a rotational spin invari-ance, among many other different phases. Recently theexistence of a different case of SSB phase (which has in-triguingly not been considered until recent years) break-ing the time translational symmetry has been under largediscussion. These phases first addressed by Wilczek in2012 [4] (and later termed as time crystals) generated anintense debate [5–10] and were soon ruled out in thermalequilibrium system by a no-go theorem [11] (for short-range interacting system) indicating in this way that theproper ground for its existence are in out-of-equilibriumconditions. In fact, theoretical studies along this direc-tion have been successful in predicting the existence oftime crystals in disparate different systems, ranging fromclosed to open systems breaking a continuous or dis-crete time translational symmetry [1, 12–32]. In particu-lar, discrete time crystal were observed experimentally in2017 in an interacting spin chain of trapped atomic ions[33] and on dipolar spin impurities in diamond [34], soonafter their theoretical preditions. Later on other sys-tems were also experimentally observed supporting suchpeculiar phases of matter [35–38]. See Refs.[39, 40] forinteresting reviews on the topics.A particularly interesting form of time crystal phasescan occur in quantum systems in contact to an environ- ment. In this case the system, also termed as bound-ary system [1], can break the time translation symmetrywhile the environment remaining time-translationally in-variant. The symmetry breaking appears only at the(macroscopic) boundary system, thus forming a so-called boundary-time crystal (BTC) (similar to surface criticalphenomena). In such phases the system shows, onlyin the thermodynamic limit, a persistent dynamics ofa macroscopic observable breaking the time translationsymmetry. A specific model supporting this phenomenol-ogy was shown in Ref.[1] composed of collective inter-acting 2-level systems. The model breaks a continuous time translational symmetry for specific regimes of inter-actions, with its collective magnetization showing peri-odic closed orbits during the dynamics.A major motivation of this work is to study differentforms BTC’s can appear in quantum systems. A betterunderstanding of its different facets and peculiar prop-erties could help us to improve our understanding forsuch non-equilibrium phases of matter and its possibleconnections to apparently different concepts [26, 41]. Inparticular, it could shed some light on the basic mecha-nisms (and limitations) supporting such phases, provid-ing insights e.g. on the role of global and dynamicalsymmetries on collective models [28], or on the nature ofits elementary excitations [43]. Moreover, we also exploreextended models with interactions between BTC phases,a subject recently investigated experimentally [44], givingrise to richer dynamical phases. We hopefully expect thata proper understanding of the fundamental properties ofsuch phases may open the way for novel applications indifferent fields, as recently proposed for the simulation ofquantum complex networks [41] or in protocols for non-Abelian braidings of Majorana edge modes in quantumcomputation [45]. a r X i v : . [ c ond - m a t . s t a t - m ec h ] F e b Figure 1: We consider collective d -level systems described by a Lindblad master equation - Eq.(1) - whose dynamics is driven bycoherent Hamiltonian terms ( ω ) and collective dissipation ( κ ). Specifically we study the cases of (A) a collective 2-level system(Eqs.(2)-(3)), (B) a pair of interacting collective 2-level systems (Eqs.(4)-(5)) and (C) a collective 3-level system (Eqs.(6)-(7)). In this work we consider different models supportingBTC’s in collective d -level systems. We first discuss indetail the collective 2-level system and further use it as abasis for the definition of the cases, featuring even richerdynamical phases. We discuss the steady state proper-ties of the collective 2-level system with a combination ofanalytical and numerical methods, characterizing theirstability from a Jacobian perspective. This allows usto make a correspondence of BTC’s to the existence of centers (also called neutrally stable fixed points) in themodel. From this perspective we trace the phase diagramof the model considering different interacting terms in theHamiltonian. We also highlight the role of the certainsymmetries in the model.In a second part of the manuscript we use these resultsas a basis for the definition of an extended model com-posed of a pair of collective interacting 2-level systems -Fig.(1)-B. We show that this model supports richer dy-namical phases. In particular we show that for certaininteractions breaking the coherent Hamiltonian symme-try (as well as apparently any quasi-conserved quantity)the model shows limit-cycles regimes, thus supporting amore robust BTC phase. Further varying the interactionthe system tends to a chaotic dynamics from subsequentperiod-doubling bifurcations. We also analyse the period-doubling bifurcation ratio of the model and its Lyapunovspectrum.In a third part of the manuscript we consider a modelwith microscopic constituents composed of 3-level sub-systems (Fig.(1)-C) rather than the 2-levels or pairs ofthem. The model deviates significantly from the sim-pler 2-level case. The collective operators now belongs toan SU (3) algebra, which cannot be reduced to SU (2) orproducts of it, which has basic implications to the globaland dynamical symmetries of the system. We study inparticular a Lindbladian with the competition of pairs of2-levels, similar to the previous case, however consideringnow the case in which they share a common level. Thephase diagram of the model shows static steady statescharacterized by a “dark level”, limit-cycle dynamics anda peculiar dynamical phase at a critical line supportingthe coexistence of multiple limit-cycles and closed orbits. This manuscript is organized as follows. In Sec.(II)we define the three different models studied in this work.We start our analysis with the simplest system, i.e., with d = 2. In Sec.(III) we derive its dynamical equationsof motion from a semiclassical approach and discuss itssymmetries and quasi-conserved quantities. We also ob-tain analytically the different steady states of the modeland analyze in Sec.(IV) their linear stability from theJacobian matrix. In Sec.(V), based on our previous re-sults, we study the effects of a specific symmetry breakingHamiltonian perturbation in the model and the stabilityof the BTC’s. In Sec.(VI) we start the study of the pairof collective systems, i.e., d = 4. We first obtain thedynamical equations of motion and discuss their symme-tries. In Sec.VII we explore the effects of interactions be-tween the pair of collective systems, showing the appear-ance of limit-cycle regimes, period-doubling bifurcationsand a route to chaos. In Sec.(IX) we move our analysisto the model with d = 3. We first introduce the basicsof SU (3) algebra and Gell-Mann basis, derive the semi-classical dynamical equations and discuss its symmetries.In Sec.(IX) we analyse the phase diagram of model. Wepresent our conclusions in Sec. (X). II. THE MODELS
In this section we define the models studied in themanuscript and its general properties. We consider mod-els with collective interactions composed of
N d -level sub-systems coupled to a Markovian environment. The timeevolution for the system is described by the master equa-tion [46], ddt ˆ ρ = ˆ L [ˆ ρ ] = i [ˆ ρ, ˆ H ] + (cid:88) i (cid:18) ˆ L i ˆ ρ ˆ L † i + 12 { ˆ L † i ˆ L i , ˆ ρ } (cid:19) , (1)with ˆ L the Lindbladian superoperator, ˆ H the coherentdriving Hamiltonian of the system and ˆ L i the Lindbladjump operators, describing the coupling of the system tothe environment. We define the three different modelsbelow, specifically, the cases with collective d = 2, 3 and4 level systems. A. Collective d = 2 -level systems We start considering the simpler case with d = 2. Inthis case the coherent Hamiltonian and Lindblad jumpoperators are defined as,ˆ H = ω ˆ S x + ω x S ( ˆ S x ) + ω z S ( ˆ S z ) , (2)ˆ L = (cid:114) κS ˆ S − , (3)where S = N/ S α = (cid:80) j ˆ σ αj / α = x, y, z are collective spin operators,ˆ S ± = ˆ S x ± iS y and ˆ σ αj are the Pauli spin operators forthe j ’th subsystem. The collective operators inherit the SU (2) algebra of their components, satisfying in this waythe commutation relations [ ˆ S α , ˆ S β ] = i(cid:15) αβγS γ ˆ S γ . Dueto the collective nature of the interactions, the modelconserves the total spin S = ( ˆ S x ) + ( ˆ S y ) + ( ˆ S z ) .The model on its simplest form, with ω x = ω z = 0, iscommonly used to describe cooperative emission in cav-ities [47–51] and was recently shown to support a timecrystal phase with the spontaneously breaking of time-translational symmetry [1]. While in the strong dissipa-tive case κ/ω > z -direction, in the weak dissipativecase κ/ω < B. Collective d = 4 -level systems In this case we consider a model describing a pair of col-lective 2-level (spin 1 /
2) systems. Specifically, we definethe coherent Hamiltonian and Lindblad jump operatorsas follows,ˆ H = ω xx S ˆ S x ˆ S x + ω zz S ˆ S z ˆ S z + (cid:88) p =1 ω x,p ˆ S xp + ω z,p ˆ S zp , (4)ˆ L = (cid:114) κ S ˆ S − , , ˆ L = (cid:114) κ S ˆ S − , , (5)where S = N/ S αp = (cid:80) j ˆ σ αj,p / p = 1 , α = x, y, z are thecollective spin operators for the p ’th collective 1 / σ αj,p are the usual Pauli spin op-erators for the j ’th spin in the p ’th collective system,and the excitation and decay operators are defined anal-ogously ˆ S ± ,p = ˆ S xp ± iS yp . The collective operators inheritthe SU (2) algebra for fixed p , while commuting other-wise: [ ˆ S αp , ˆ S βp (cid:48) ] = iδ p,p (cid:48) (cid:15) αβγS γ ˆ S γp . Due to the collectivenature of the model, it conserves the total spin for each collective spin system S p = ( ˆ S xp ) + ( ˆ S yp ) + ( ˆ S zp ) for p = 1 , ω xx = ω zz = 0 there is no coupling be-tween the two collective systems and the physics reducesto the simpler d = 2 case. On the other hand, if the cou-pling between the collective systems in nonzero, as e.g. ω xx (cid:54) = 0, one may expect the strengthen of the persistentoscillations in the time crystal phase, since such couplingscan induce local spin excitations on each collective systemthus enhancing the effect of coherent collective drivingsˆ S x , . We will discuss in more detail, in the next sections,the effects of the different terms in the model leading toa richer phase diagram. C. Collective d = 3 -level systems In this case the model describe cooperative evolutionof a collection of three-level subsystems ( d = 3). Westudy how a pair of collective 2-level subsystems compete,or hybridize, when they share a common energy level.Specifically, we study the competition of two dissipativechannels with the Lindbladian given as follows,ˆ L = (1 − δ ) ˆ L + δ ˆ L , (6)where 0 ≤ δ ≤ L mn acts only in the pair of lev-els m and n - see Fig.(1). The Lindbladians ˆ L mn are de-fined similarly to the d = 2 case, with coherent Hamilto-nian ( ˆ H ( mn ) ) and Lindblad jump operator ( ˆ L ( mn ) ) givenby, ˆ H ( mn ) = ω mn ˆ S xmn , ˆ L mn = (cid:114) κ mn S ˆ S − ,mn . (7)where S = N/ S αmn = (cid:80) Nj =1 ˆ σ αj,mn / α = x, y, z and m, n = 1 , , m, n ) lev-els. The operators ˆ σ αj,mn are the usual Pauli spin oper-ators for the j ’th subsystem in the pair of ( m, n ) levels.The collective excitation and decay operators are definedanalogously, ˆ S ± ,mn = ˆ S xmn ± iS ymn .Similar to the previous d = 4 case, the model here con-siders the competition of a pair of collective two-level sub-systems. A major contrast comes however from the factthat, due to the shared collective level in the d = 3 case,the collective operators form an SU(3) algebra, whichcannot be reduced to an SU(2) as in the d = 2 case,neither to a pair SU(2) ⊗ SU(2) as in the d = 4 case.The dynamics are thus expected to be different from theprevious cases, e.g one can already notice that the totalspin for each pair of two levels is not conserved anymore.The conserved quantities in this case are rather differ-ent, given by the two independent Casimir operators ofthe algebra, a quadratic and cubic operator, respectively.While the quadratic can be seen as a vector norm in thespace of group operators, the cubic operator is rathernon intuitive. We discuss in more detail these operatorsin Sec.(IX). III. d = 2 : DYNAMICAL EQUATIONS OFMOTION, SYMMETRIES AND STEADY STATES We study in this section the dynamical equations ofmotions for the collective 2-level system, their symme-tries and (quasi)-conserved quantities, as well as obtainthe steady states of the model for varying couplings inthe Lindbladian.
Dynamical Equations of Motion.-
In order to studythe dynamics and steady state properties of the model wederive the dynamical equations of motion for its collectiveobservables and employ a semiclassical approach, usualto collective systems. The dynamics of a general operatorˆ O within the Heisenberg picture is given by, d (cid:104) ˆ O (cid:105) dt = i (cid:104) [ ˆ H, ˆ O ] (cid:105) + (cid:88) i (cid:104) [ ˆ L † i , ˆ O ] ˆ L i + ˆ L † i [ ˆ O, ˆ L i ] (cid:105) . (8)Considering the collective 1 / S α andusing their SU(2) commutation relations we find the cor-responding dynamical equations, ddt (cid:104) ˆ S x (cid:105) = − ω z S (cid:16) (cid:104) ˆ S z ˆ S y (cid:105) + (cid:104) ˆ S y ˆ S z (cid:105) (cid:17) + κ S (cid:16) (cid:104) ˆ S z ˆ S x (cid:105) + (cid:104) ˆ S x ˆ S z (cid:105) + (cid:104) ˆ S x (cid:105) (cid:17) ,ddt (cid:104) ˆ S y (cid:105) = − ω (cid:104) ˆ S z (cid:105) + ω z − ω x S (cid:104) ˆ S x ˆ S z + ˆ S z ˆ S x (cid:105) + κ S (cid:16) (cid:104) ˆ S z ˆ S y (cid:105) + (cid:104) ˆ S y ˆ S z (cid:105) − (cid:104) ˆ S y (cid:105) (cid:17) , (9) ddt (cid:104) ˆ S z (cid:105) = ω (cid:104) ˆ S y (cid:105) + ω x S (cid:16) (cid:104) ˆ S y ˆ S x (cid:105) + (cid:104) S x ˆ S y (cid:105) (cid:17) − κS (cid:16) (cid:104) ˆ( S y ) (cid:105) + (cid:104) ( ˆ S x ) (cid:105) + (cid:104) ˆ S z (cid:105) (cid:17) . It is convenient to define the reduced operator ˆ m α =ˆ S α /N . These operators commute in the thermodynamiclimit [ ˆ m α , ˆ m β ] = i(cid:15) αβγ ˆ m γ /N motivating us to performa second order cumulant approximation (usually calledas semiclassical approach) for their expectation values: (cid:104) ˆ m α ˆ m β (cid:105) ∼ = (cid:104) ˆ m α (cid:105)(cid:104) ˆ m β (cid:105) . In this way in the thermodynamiclimit the dynamical equations of motion are closed andgiven by the following set of nonlinear differential equa-tions: ddt m x = m z ( − ω z m y + κm x ) ,ddt m y = m z (2( ω z − ω x ) m x − ω + κm y ) , (10) ddt m z = ω m y − κ (( m x ) + ( m y ) ) + 2 ω x m x m y . where we use m α ≡ (cid:104) ˆ m α (cid:105) to simplify our notation. Symmetries and (Quasi) Conserved Quantities.-
Wefirst see that these dynamical equations conserves the to-tal spin of the system N = ( m x ) + ( m y ) + ( m z ) -as expected since it should accurately describe the Lind-bladian dynamics which has explicitly such a symmetry. Figure 2: Steady states of the model with m z = 0: Thecolored curves are the algebraic condition of Eq.(14) whilethe black circle the normalization condition ( m x ) + ( m y ) +( m z ) = 1. The physical steady states correspond to theintersection of the two curves. We show different cases ofsystem parameters, highlighting the cases without such steadystates, a pair and two pairs. Moreover, the dynamical equations also have a quasi-conserved quantity for ω z > ω x given by, R = ( − iκ + 2 ω z ) log (cid:18) im x + m y − ω ( κ − iω z ) (cid:19) − ( − iκ − ω z ) log (cid:18) − im x + m y − ω ( κ + 2 iω z ) (cid:19) . (11)Due to the logarithmic function this quantity is definedup to integer multiples of 2 κπ , and that is the reason weprefer to define it as a quasi-conserved quantity.It is worth noticing that the dynamical equations havea reversibility symmetry, given by the following transfor-mation t → − t, m x → m x , m y → m y , m z → − m z . (12)Furthermore we see a specific structure for the m x and m y dynamical equations, where the term m z can be fac-tored out leading to specific conditions for a steady state. Steady States.-
We focus now on the analysis of thesteady states of the semiclassical dynamical equations,and obtain results both analytical as algebraic (quasi-analytical), depending on the fixed points. In order toobtain the steady states of the system we must solve dm α /dt = 0 for α = x, y, z . Recalling the specific struc-ture for the dynamical equations of motion, with the m z term factoring out, we can consider two different casesfor the steady states: (i) m z (cid:54) = 0 and (ii) m z = 0. Let usstart analyzing the first case.(i) Considering m z (cid:54) = 0 the solutions for the dynamical Figure 3: Azimuthal angle for the steady states with m z = 0(Eq.(14)), for varying ω x field and system parameters (a) κ = 0 . , ω = 1, and (b) κ = 1 , ω = 0 .
5. The steady statesare independent of the ω z field. equations correspond to the pair of fixed points, x ∗ = 2 ω z ω ω z ( ω z − ω x ) + κ ,y ∗ = κω ω z ( ω z − ω x ) + κ , (13) z ∗ = ± (cid:115) − ω ω z + κ (4 ω z ( ω z − ω x ) + κ ) , where we denote as ( x ∗ , y ∗ , z ∗ ) the fixed points for sim-plicity of notation. In order for these fixed points bea physical steady state they must also satisfy the normcondition of fixed total spin N . The system can supportin this way zero or one pair of steady states with m z (cid:54) = 0.(ii) Considering the case of m z = 0 the solutions forthe dynamical equations correspond now to the algebraiccondition, y ∗ = κω + 2 ω x x ∗ , (14)and we recall that this fixed point solutions are physi-cal steady state only if this geometrical place intersectsthe fixed norm circle ( m x ) + ( m y ) = N , thus satisfy-ing the conservation of the total spin. In this case thesteady states come in pairs, as in the previous case, how-ever we can have 0 , κ, ω and ω x - see Fig.(2) for illustrativecases. We also show in Fig.(3) the steady states mag-netization for varying Hamiltonian parameters. Since m z = 0 we use the azimutal angle in the x − y plane, φ = arctan ( y ∗ /x ∗ ), in order to completely describe thesteady state. The steady state magnetization is indepen-dent of the ω z field (Eq.(14)). We see that varying the ω x field, both starting from a trivial ( κ/ω >
1) or time-crystal phase ( κ/ω < x − y plane. IV. STABILITY ANALYSIS
In this section we perform a stability analysis for thesteady states of the model. A simple approach is basedon the linearization of dynamical equations of motionaround the fixed point, which is effectively described by the Jacobian matrix. The spectrum of the Jacobian pro-vides information on the stability of the fixed points.Specifically, the dynamical equations of motion of thesystem can be written as dm α /dt = f α ( m x , m y , m z ),with α = x, y, z and f α a nonlinear function on thevariables. We can define the displaced variables u α = m α − α ∗ around the fixed points and perform a seriesexpansion. We obtain that, (cid:18) du x dt , du y dt , du z dt (cid:19) T = ˆ J ( u x , u y , u z ) T + O (( u α ) , u α u β ) , (15)where (cid:126)v T denotes the transposed of the line vector (cid:126)v ,ˆ J is a 3 × J ) αβ = ∂f α /∂β and the correction terms O (( u α ) , u α u β ) are quadratic on the displaced variablesand can be neglected within a linear approximation. Theeffective dynamical equations of motion around the fixedpoint are in this way linear differential equations whichcan be solved by the eigenspectrum of the Jacobian ma-trix.In case the eigenvalue of the Jacobian matrix have neg-ative (positive) real part the fixed point is an attractor(repeller) and robust to nonlinear terms in the above se-ries expansion. In case where the eigenvalues have also animaginary term the dynamics have spirals towards (away)to the fixed point. Fixed points with nonnull eigenvaluereal part are usually called as hyperbolic fixed points.Another interesting case is when the eigenvalue of theJacobian matrix is purely imaginary. In this case thefixed point is denoted as a center . We will see that theexistence of centers is intimately related to the presenceof time crystal behavior in this model. We recall thatthe robustness of centers to nonlinear terms is in princi-ple not guaranteed, as will become clear when studyingspecific symmetry breaking Hamiltonian perturbations inSec.(V).We start our analysis computing the Jacobian matrixfor the fixed points of the model, obtaining that,ˆ J = J J J J J J J J J , (16)where, J = J = κz ∗ , J = − ω z z ∗ ,J = − ω z y ∗ + κx ∗ , J = 2( ω z − ω x ) z ∗ , (17) J = 2( ω z − ω x ) x ∗ − ω + κy ∗ , J = − κx ∗ + 2 ω x y ∗ ,J = − κy ∗ + 2 ω x x ∗ + ω , J = 0 . Using sympy package from python we obtain that the de-terminant of this matrix for the pair of fixed points with z ∗ (cid:54) = 0 (Eq.(13)) is zero. The same is true for any fixedpoint with z ∗ = 0 in which case the matrix is stronglysimplified having the first 2 × Case z ∗ (cid:54) = 0 . For the pair of fixed points z ∗± of Eq.(13) we computethe eigenvalues using sympy package from python. Weobtain that their eigenvalues are given by, λ [ z ∗± ] = (cid:26) , z ∗± ± √ ω z Aκ − ω x ω z + 4 ω z (cid:27) , (18)where, A = (cid:0) κ ω x + κ ω ω z + 16 κ ω x ω z + 4 ω ω z + 16 ω x ω z +48 ω x ω z (cid:1) − (cid:0) κ ω z + κ ω ω x + 8 κ ω z + 8 κ ω x ω z +4 ω ω x ω z + 48 ω x ω z + 16 ω z (cid:1) . (19)Apart from the trivial null eigenvalue, we see that thenontrivial eigenvalue has always a real part. As a con-sequence the possible steady states of the model with z ∗ (cid:54) = 0 are always hyperbolic fixed points. Case z ∗ = 0 . For the case with z ∗ = 0 the Jacobian matrix is sim-plified and one can obtain analytically their eigenvalues,which have the form, λ = (cid:110) , ±√ B (cid:111) , (20)where, B = (cid:0) κω y ∗ + 8 κω x x ∗ y ∗ + 2 ω ω z x ∗ + 4 ω x ω z ( x ∗ ) (cid:1) − (cid:0) ω + 2 κ N + 4 ω ω x x ∗ + 4 ω x ( x ∗ ) + 4 ω x ω z ( y ∗ ) (cid:1) . (21)In this case we see that the eigenvalues (apart from thetrivial one) are purely real or imaginary, thus corre-sponding either to hyperbolic fixed points or centers. Phase Diagram From Stability Analysis.-
We studythe full phase diagram of the model from the perspec-tive of their Jacobian stability. In the simpler case of themodel, with ω x = ω z = 0, we know the system supportsa trivial phase in the strong dissipative case ( κ/ω > κ/ω < λ ± (Eqs.(18), (20))and steady states ( x ∗ , y ∗ , z ∗ ) (Eqs(13) and (14)): • Strong dissipative : λ ± = ± (cid:112) κ − ω , (cid:18) , ω κ , ± (cid:113) − (cid:0) ω κ (cid:1) (cid:19) . (22) • Weak dissipative :: λ ± = ± (cid:112) κ − ω , (cid:32) ± (cid:114) − (cid:16) κω (cid:17) , κω , (cid:33) . (23)While in the strong dissipative case one has purelyreal Jacobian eigenvalues, thus characterizing hyperbolicsteady states with a pair of an attractor and a repeller,in the weak dissipative case we have purely imaginaryeigenvalues, corresponding to centers. We can thus asso-ciate the presence of centers in the model to the existenceof closed orbits and a time crystal behavior.From this perspective we map the possible existenceof time crystals in the model to the presence of steadystate centers. We analyze the Jacobian spectrum for thefull phase diagram (varying ω x and ω z couplings) consid-ering both strong and weak dissipative cases. We showour results in Fig.(4). We see that while small couplingsmay destabilize (or shrink the phase space support for)the time crystal phase, strong couplings can stabilize itboth in the weak as in the strong dissipative case, lead-ing even to regions with the presence of four stable cen-ters in the model. We show in Fig.(5) phase portraits( Q, P ) defined by m z = Q , m x = (cid:112) − Q cos(2 P ) and m y = (cid:112) − Q sin(2 P ) for a few points of the phase di-agram, corroborating our association for the presence ofcenters with time crystals in the model. V. SYMMETRY BREAKING PERTURBATION
In this section we study the effects of a specific pertur-bation in the system breaking the coherent Hamiltonian Z symmetry. Motivated by the results of Sec.(IV) show-ing that centers can only occur for the steady states with m z = 0, we consider as perturbation a field along the z - direction. In this way the model Hamiltonian is givenby ˆ H δ z = ˆ H + δ z ˆ S z and the corresponding semiclassicaldynamical equations are obtained:( ddt m x ) δ = ( ddt m x ) δ =0 − δm y , (24)( ddt m y ) δ = ( ddt m y ) δ =0 + δm x , (25)( ddt m z ) δ = ( ddt m z ) δ =0 , (26)where ( dm α /dt ) δ =0 denotes the unperturbed semiclassi-cal equations (Eq.(11)). The Jacobian of the system canbe written as ˆ J δ = ˆ J + δ ˆ V , (27)
Figure 4:
Phase diagram from Jacobian analysis: we showthe different phases of the model from the perspective of theirsteady states stability, for varying Hamiltonian couplings ω z and ω z . We show in (a) our results for a system in the weakdissipative case, with fixed κ = 0 . , ω = 1 and in (b) forthe strong dissipative case, with fixed κ = 1 , ω = 0 .
5. Thesteady states of the model can be centers or single hyperbolicsteady states with z ∗ = 0, or pairs of attractor-repulsor hyper-bolic steady states with z ∗ (cid:54) = 0. We thus describe the phasesby the 3-vector ( n c , n sh , n ph ) where n c denotes the numberof such centers, n sh the number of single hyperbolic steadystates and n ph the number of attractor-repulsor pairs. Thered lines highlight when a pair of attractor-repulsor are cre-ated/annihilated, the vertical black lines correspond to thecreation/annihilation of a pair of steady states with z ∗ = 0 -which can come as centers or single hyperbolic steady states- and the dark curves represent a change of stability of centerto single hyperbolic steady states ( n c ↔ n sh ). where ˆ J is the unperturbed Jacobian (Eq.(16)) andˆ V = − , (28) is the perturbation matrix. We can study the effects ofsuch perturbation in the spectral properties of the Jaco-bian using Perturbation Theory for general matrices [53](notice here that the Jacobian matrix is not an Hermi-tian matrix). The idea follows similarly to the simplerHermitian case. Any matrix has generalized right ( (cid:126)u i )and left ( (cid:126)w i ) eigenvectors defined as,ˆ J(cid:126)u i = λ i (cid:126)u i , ˆ J † (cid:126)w i = λ ∗ i (cid:126)w i , (29)where left and right eigenvectors obey the ortogonal-ity property Tr( (cid:126)u † i (cid:126)w j ) ∝ δ ij , and λ i are the general-ized eigenvectors. We assume that the eigenvalues arenon-degenerated, and thus we are dealing with non-degenerate Perturbation Theory. Expanding the eigen-values and eigenvectors in terms of the perturbative term δ , we can write then as, λ i = ∞ (cid:88) j =0 δ j λ ( j ) i , (cid:126)u i = ∞ (cid:88) j =0 δ j (cid:126)u ( j ) i , (cid:126)w i = ∞ (cid:88) j =0 δ j (cid:126)w ( j ) i , (30)where the index j denotes the correction order in pertur-bation theory, i.e. j = 0 corresponds to the eigenvaluesand eigenvectors of the unperturbed Jacobian J . Com-bining Eq.(27) and Eq.(30) into Eq.(29), one can findthe recursive expressions for general j ’th order correc-tions. In particular, the corrections to the eigenvalues ofthe Jacobian are given by, λ ( j ) i = Tr (cid:16) (cid:126)w (0) † i ˆ V (cid:126)u ( j − i (cid:17) − j − (cid:88) k =1 λ ( k ) i Tr (cid:16) (cid:126)w (0) i (cid:126)u ( j − k ) i (cid:17) . (31)Focusing on the first order terms ( j = 1) we find thatthe perturbative eigenvalue corrections λ (1) i = Tr( (cid:126)w (0) † i ˆ V (cid:126)u (0) i ) (32)for the steady states with z ∗ = 0 (Eq.(14)) in our modelare given by, λ (1)1 = − κx ∗ − ω x y ∗ )( − κy ∗ + ω + 2 x ∗ ( ω x − ω z )) + ( κx ∗ − ω z y ∗ )( − κy ∗ + ω + 2 ω x x ∗ )2( κx − ω x y ∗ )( κx ∗ − ω z y ∗ ) ,λ (1)2 = 2( κx ∗ − ω x y ∗ )( − κy ∗ + ω + 2 x ∗ ( ω x − ω z )) + ( κx ∗ − ω z y ∗ )( − κy ∗ + ω + 2 ω x x ∗ )2 κ x ∗ + 2 κ y ∗ − κω y ∗ − κω x x ∗ y ∗ + ω + 4 ω ω x x ∗ − ω ω z x ∗ + 4 ω x − ω x ω z x ∗ + 4 ω x ω z y ∗ , (33) λ (1)3 = 2( κx ∗ − ω x y ∗ )( − κy ∗ + ω + 2 x ∗ ( ω x − ω z )) + ( κx ∗ − ω z y ∗ )( − κy ∗ + ω + 2 ω x x ∗ )2 κ x ∗ + 2 κ y ∗ − κω y ∗ − κω x x ∗ y ∗ + ω + 4 ω ω x x ∗ − ω ω z x ∗ + 4 ω x − ω x ω z x ∗ + 4 ω x ω z y ∗ . Figure 5: Phase space portraits (
Q, P ) for the collective 2-level system with κ = 0 . , ω = 1 and different Hamiltonian couplings ω x , ω z . We show in (a-b) the results for a fixed ω x = 0 with (a) ω z = 0 and (b) ω z = 2. In (c-f) we show the phase portraitsfor the phase diagram line with fixed ω x = 2, where (c) ω z = − (d) ω z = 0, (e) ω z = 1, (f ) ω z = 3. Since x ∗ , y ∗ ∈ (cid:60) the first order corrections are purelyreal terms, implying that the steady state centers be-come hyperbolic steady states. In this way, the closedorbits characteristics of time crystals are destroyed andwe have instead spirals towards or away from the fixedpoints of the model, with characteristic times capturedby the real eigenvalues λ (1) i . These results follow in accor-dance with closely related p, q - interacting model recentlystudied in Ref.[27, 28], shown also to support BTC onlyin the absence of a Z Hamiltonian symmetry breakingperturbation.
VI. d = 4 : DYNAMICAL EQUATIONS OFMOTION AND SYMMETRIES We move our studies now to the case of a pair of col-lective 2-level (1 / m αp = ˆ S αp /N and closing the expectations values in thesecond cumulant (cid:104) ˆ m αp ˆ m βp (cid:105) ∼ = (cid:104) ˆ m αp (cid:105)(cid:104) ˆ m βp (cid:105) we obtain the semiclassical dynamical equations of motion: ddt m x = − ω zz m y m z − ω z, m y + κ m x m z ,ddt m y = ω z, m x − ω x, m z − ω xx m z m x + ω zz m x m z + κ m y m z ,ddt m z = ω x, m y + ω xx m y m x − κ (cid:0) ( m x ) + ( m y ) (cid:1) ,ddt m x = − ω zz m z m y − ω z, m y + κ m x m z , (34) ddt m y = ω z, m x − ω x, m z − ω xx m x m z + ω zz m z m x + κ m y m z ,ddt m z = ω x, m y + ω xx m x m y − κ (cid:0) ( m x ) + ( m y ) (cid:1) . Symmetries and Conserved Quantities.-
The dynam-ical equations conserve the total spin for each collective1 / N p = ( m x ) p +( m y ) p +( m z ) p for p = 1 , ω xx , ω zz do not breakthe reversibility symmetry of the equations. In particu-lar, in the case where ω z, = ω zz = 0 the equations stillhave the factorization structure for the m z terms in the m x and m y dynamical equations. In this case one canproceed the analysis similarly to Ref.[1] and show thatthe system do have (quasi) conserved dynamical quanti-ties. In case ω z, or ω zz (cid:54) = 0, however, the coupling de-stroys this simpler factorization structure making incon-clusive the existence of conserved quantities (one breaksalso the Z Hamiltonian symmetry). One could consider,however, the simpler case with only ω z, = 0 and equallocal couplings for both collective spins ( ω x, = ω x, , κ = κ ) and study the specific case where the collec-tive spins are initially the same in the evolution. In this t m i z m z m z t m i z m z m z t m i z m z m z
100 105 110 115 120 125 130 t m i z m z m z Figure 6: We show the dynamics of the spin magnetizationalong the z-axis for system parameters of Eq.(36) and dif-ferent couplings ω zz with (cid:126)m (0) = (0 , , , , , (upper-left panel) a steady statewith fixed magnetization for ω zz = 0, (upper-right panel) aBTC with limit-cycles for ω zz = 0 . (bottom-left panel) achaotic dynamics for ω zz = 0 .
58 and (bottom-right panel) limit-cycles with period doubling oscillations for ω zz = 1. case they shall also remain the same throughout all thedynamics, ˆ m α ( t ) = ˆ m α ( t ) ∀ t , and we recover the fac-torization structure in the dynamical equations, thus theexistence of (quasi) conserved quantities. VII. ROUTE TO CHAOS
In the general case of nonzero couplings in the modelthe absence of (quasi) conserved quantities - beyond thenorm N of the collective spins - and no Z Hamiltoniansymmetry turns the analysis of the steady states and dy-namics of the system more intricate. On the other hand,it allows the possibility of more complex dynamics withricher dynamical phases. A particularly interesting caseoccurs when, (i) w x,p /κ p > , (ii) w z,p /κ p (cid:28) , (iii) ω xx /κ p (cid:29) , (35)(iv) ω zz /κ p (cid:54) = 0 , for p = 1 ,
2. Conditions (i) and (ii) cannot alone sup-port BTC’s, as shown in the previous sections, ratherthey are characterized by hyperbolic steady states with z ∗ (cid:54) = 0. Condition (iv) can correlate these steady stateswith (iii) inducing coupled spin excitations on the system.The system in this case has no (quasi) conserved quan-tities, and the appearance of BTC’s shall be due to the collective dynamics of both spin systems ( hybridization ).Specifically, we study the case with system parametersgiven by (i) w x,p = 2 , p = 1 , , (ii) w z, = 0 . , w z, = 0 . , (iii) w xx = 3 , (36)with κ = κ = 1 and for varying coupling ω zz . We showin Fig.(6) the dynamics for the magnetization along the z -axis for different cases of the coupling ω zz . We see thatthe system can support different phases ranging from (i)a ferromagnetic phase with nonzero magnetization for itssteady states, (cid:104) ˆ m zp ( t → ∞ ) (cid:105) (cid:54) = 0, (ii) BTC’s characterizedby limit-cycle oscillations, where after an initial transienttime the magnetization oscillates in a given orbit indefi-nitely in time; moreover, further increasing the couplingthe limit-cycles oscillations are followed by period dou-bling bifurcations till (iii) reaching a chaotic dynamics.Interestingly, we also see (iv) 3-cycle periodic windowsintercalated by the stable period doubling bifurcationsand the chaotic regime.A global picture of the dynamical phases in the modelare shown in the orbit diagram of Fig.(7). The orbit di-agram corresponds to the local minimums in the timeevolution of the magnetization along the z-axis, obtainedafter an initial transient time. The transient time is re-lated to the relaxation of the initial collective spin statetowards the limit-cycle orbits, static steady states orchaotic regime. In our numerical simulations for the orbitdiagram we used a fixed initial state for the dynamics,with m y (0) = m z (0) = 1 and zero otherwise. We ob-served however that the orbit diagram is qualitativelysimilar considering a few different initial conditions [54].We can accurately determine in the orbit diagram thefirst period doubling bifurcations as we increase the cou-pling ω zz . We obtain a bifurcation ratio b ≈ .
2, where b n = ω ( n − zz − ω ( n − zz ω ( n ) zz − ω ( n − zz , (37)with ω ( n ) zz the coupling corresponding to the n ’th perioddoubling bifurcation in the orbit diagram. Interesting tocompare with Feigenbaum constant for the seminal logis-tic map, in which one has b n →∞ ≈ .
67. In our modelwe obtained a slightly different value, which could indi-cate a different universality class for the period doublingcascades towards chaoticity. We remark, however, thatwe were able to obtain only the first bifurcation ratio b (not precisely n → ∞ ) so an irrefutable conclusion onthe universality class cannot be drawn at this moment.It remains as a very interesting perspective for a futurework.We also study the Lyapunov exponents in the limit-cycle and chaotic regimes [52, 55–57]. We obtain the fullLyapunov spectrum of the dynamics, describing the meangrowth of an n -dimensional volume ( n = 6 in our case)in the tangent space. We use Benettin’s approach in our0 Figure 7: (top panel)
Orbit diagram for system parametersof Eq.(36) and varying ω zz . We show in the bottom panel a zoom for the region around ω zz ≈ .
37. The orbit diagramis obtained from the local minimums in the time evolution ofthe system magnetization along the z-axis. We see a multi-tude of dynamical phases. For ω zz (cid:46) . . (cid:46) ω zz (cid:46) . ω zz ∼ .
53 we see the appearance of a 3-cycleperiodic window. For larger ω zz (cid:38) analysis, i.e., employing recursively the Gram-Schmidtorthonormalization procedure during the stretching andfolding of the n -dimensional volume. In our numericalsimulations we set the same initial state as in the or-bit diagram, and perform an initial transient evolutionup to t = 500. After this initial evolution the collectivespin is close to their corresponding dynamical phases, andwe employ Bennetin’s algorithm to obtain the Lyapunovspectrum. We show our results in Fig.(8). We see thatwhile for ω zz = 0 .
35 the spectrum is all non-positive, inthe chaotic region with ω zz = 0 .
45 the largest Lyapunovexponent is positive, corroborating our orbit diagram ex-pectations of a limit-cycle and chaotic regimes, respec-tively. We further notice that the sum of the Lyapunovspectrum is negative in both cases, typical of dissipativedynamical equations. t -1-0.500.51 Λ L y a p j t -1-0.500.51 Λ L y a p j Figure 8: We show the Lyapunov spectrum Λ lyap j obtainedfrom Benettin’s approach. We see the convergence of theΛ lyap j exponents towards its assymptotic value ( t → ∞ ). Wealso see the convergence of the assymptotic exponents as wedecrease the time steps. In particular we obtain that thelargest Lyapunov exponent for ω zz = 0 .
35 approaches zeroin the long time limit as a power law with time, while for ω zz = 0 .
45 it approaches approximately 0 .
15. We used a timestep dt = 10 − in our numerical simulations. VIII. d = 3 : DYNAMICAL EQUATIONS OFMOTION AND SYMMETRIES In this section we analyse the model with collective d = 3-level system. In this case it is convenient to workwithin the Gell-Man basis for the three-level subsystems,ˆ g j = , ˆ g j = − i i , ˆ g j = − , ˆ g j = , ˆ g = − i i , ˆ g j = , ˆ g j = − i i , ˆ g j = √ − (38)corresponding to an Hermitian basis for the j ’th sub-system. The collective operators ˆ G k = (cid:80) Nj =1 ˆ g kj inheritdirectly the algebra of their microscopic constituents, i.e. the SU(3) algebra of the Gell-Man basis, given by,[ ˆ G a , ˆ G b ] = i (cid:88) c =1 f abc ˆ G c , (39)with a, b = 1 , ..., f abc the structure constant totallyantisymmetric under the exchange of any pair of indices.The explicit form of the non-zero SU (3) structure con-stants are listed in the Table I. Any collective operator1 Table I: Non-zero structure constants f abc of SU (3) abc f abc abc f abc
123 1 345 − − √ √ Table II: Non-zero structure constants d abc of SU (3) abc d abc abc d abc √ − − √ − √ − − √ − √ √ − √ − √ can be decomposed in this basis. Specifically, the coher-ent Hamiltonian terms of the model are decomposed asˆ S x = ˆ G , ˆ S x = ˆ G , while the decay operators are givenby ˆ S ± , = ˆ G ± i ˆ G and ˆ S ± , = ˆ G ± i ˆ G .We can also define number operators for the three col-lective energy levels, corresponding to their collective oc-cupation, as follows,ˆ N = 13 I + ˆ G + 1 √ G , ˆ N = 13 I − ˆ G + 1 √ G , (40)ˆ N = 13 I − √ G , where I is the identity operator. Dynamical equations of motion.
Following the sameapproach as in the previous sections, we define the oper-ators ˆ m k = ˆ G k / ( N/
2) and close the expectations valuesin the second cumulant (cid:104) ˆ m k ˆ m (cid:96) (cid:105) ∼ = (cid:104) ˆ m k (cid:105)(cid:104) ˆ m (cid:96) (cid:105) . We ob-tain the following semiclassical dynamical equations ofmotion in the thermodynamic limit ( N → ∞ ): dm dt = k (1 − δ ) m m + 12 δ (cid:0) ω m − k (cid:0) m m + m m (cid:1)(cid:1) , (41) dm dt = (1 − δ ) (cid:0) − ω m + k m m (cid:1) + 12 δ (cid:0) − ω m + k (cid:0) m m − m m (cid:1)(cid:1) , (42) dm dt = (1 − δ ) (cid:0) ω m − k (cid:0) ( m ) + ( m ) (cid:1)(cid:1) + 12 δ (cid:0) − ω m + k (cid:0) ( m ) + ( m ) (cid:1)(cid:1) , (43) dm dt = −
12 (1 − δ ) (cid:0) ω m + k (cid:0) m m − m m (cid:1)(cid:1) + 12 δ (cid:0) ω m + k (cid:0) m m − m m (cid:1)(cid:1) , (44) dm dt = 12 (1 − δ ) (cid:0) ω m − k (cid:0) m m + m m (cid:1)(cid:1) + 12 δ (cid:0) − ω m + k (cid:0) m m + m m (cid:1)(cid:1) , (45) dm dt = 12 (1 − δ ) (cid:0) − ω m + k (cid:0) m m + m m (cid:1)(cid:1) + 12 δk (cid:16) √ m m − m m (cid:17) , (46) dm dt = 12 (1 − δ ) (cid:0) ω m + k (cid:0) m m − m m (cid:1)(cid:1) + 12 δ (cid:16) ω (cid:16) m − √ m (cid:17) + k (cid:16) √ m m − m m (cid:17)(cid:17) , (47) dm dt = √ δ (cid:0) ω m − k (cid:0) ( m ) + ( m ) (cid:1)(cid:1) . (48) Symmetries and conserved quantities.
Since all collec-tive operators can be decomposed in the SU (3) basis, themodel have conserved quantities given by the Casimir el-ements of the algebra (the element which commutes withall operators of the group). The two independent Casimirelements in the SU (3) algebra correspond to a quadratic ( ˆ C ) and a cubic operator ( ˆ C ), defined asˆ C = (cid:88) j =1 (cid:0) ˆ m j (cid:1) , (49)ˆ C = (cid:88) a,b,c d abc ˆ m a ˆ m b ˆ m c , (50)respectively, where d abc is symmetric under the inter-change of any pair of indices. The explicit form of the2 δ ω / κ t n i ( t ) t n i ( t ) t n i ( t ) ω/ κ = 0 . , δ = 0 . ω/ κ = 0 . , δ = 0 . ω/ κ = 0 . , δ = 0 . CriticalLine
Phase I Phase IIPhase III(BTC)
Figure 9: Phase diagram for the collective 3-level system - Eqs.(6),(7) with couplings of Eq.(51). We show three differentphases supported by the model, phase I, II and III, with insets as illustrative of their dynamics. The initial state for thedynamics is given by the collective occupation of a single level n ( t = 0) = 1. While phases I and II show static steadystates with different characteristics, phase III features BTC’s with limit-cycle dynamics. We highlight in red the critical lineat δ = 1 /
2, supporting multiple limit-cycle attractors as well as closed orbit dynamics, as we discuss in the main text.. non-zero d abc are listed in Table II. We see that C can beseen as a conservation of the norm in the basis of collec-tive operators. This is similar to the d = 2-dimensionalcase where the norm is associated to the total spin ofthe system, represented as the surface of a 3-dimensionalsphere. In this case, however, the basis has a higher di-mensionality, and the norm is related to the surface ofa 8-dimensional hypersphere. The interpretation of thissurface with the the total spin of the system is not di-rect anymore. This is also the case of the second Casimirelement, C . It is cubic operator in the operator basis,further constraining the dynamics on the surface of the8-dimensional hypersphere. An exact interpretation ofsuch constraint is, however, also not clear on physicalgrounds.Moreover, even though each pair of energy levelsin the subsystems are coupled similar to the d = 2-dimensional case, the symmetries and conserved quan-tities present there are no longer present in this case,namely: reversibility (Eq.(12)), quasi-conserved quanti-ties R (Eq.(11)) and total spin N . IX. d = 3 : PHASE DIAGRAM We study in this section the phase diagram for the d = 3 model. We focus in the case of the Lindbladianwith couplings ω = αω ≡ ω, κ = ακ ≡ κ, (51)with α ∈ (cid:60) . We set α = 1 in all of our analysis forsimplicity, since different values corresponds simply to a renormalization of the Lindbladian and consequently its δ parameter. Specifically, for a different α ’ coupling wesee that L α,δ = c L α (cid:48) ,δ (cid:48) where c = δδ (cid:48) and δ (cid:48) = δα (cid:48) / ( α (1 − δ ) + δα (cid:48) ), thus both Lindbladians share the same steadystates and attractors of the dynamics.The phase diagram is shown in Fig.(9), characterizedby different static as well as time crystal phases. In allof our analysis we perform the dynamics starting froma few different initial state conditions and study its evo-lution towards their corresponding dynamical or staticattractors. Precisely, we consider initial states with thecollective system occupying the same level, i.e., with oc-cupation number (cid:104) n i ( t = 0) (cid:105) = 1 for i = 1 , δ = 1 / δ = 0 or 1,respectively, as we will discuss in more detail. We remarkhowever that since we deal with an 8-dimensional space,we do not preclude the existence of different attractorsin the model. A throughout analysis of the full set ofsteady states and their dependence on the initial condi-tions remains as an interesting perspective. We discussin detail the different phases of model below. • Extremal lines: for couplings δ = 0 or 1 one recoversthe SU (2) model of Sec.(III) for the pair of energy levels(1 ,
2) or (2 , ω/κ . • Phase I: for couplings 0 < δ < /
2, the Lindbaldianacts stronger on levels (1 , L , , even if small, tends to destroythe SU (2) organization on these levels. We obtain the3 Figure 10: We show the occupation numbers for the three col-lective energy levels of the system: (left panel) steady statesof Phase I - Eq.(52) - for 0 < δ < / , ∀ ω/κ ; (right panel)steady states of Phase II - Eq.(53) - for 1 / < δ < ω/κ < / following steady state attractor for the model, m = √ C δ δ − δ + 1 m = − − δ ) δ m m = 1 √ − − δ ) δ ) m m i = 0 , i = 1 , , , , . (52)where C = 4 / m i = 0 , i = 1 , , , ω/κ and the phasetransition at δ = 1 / n and n equilibrate.The conditions of Eq.(52) is actually valid as a steadystate ∀ δ , i.e., it is a solution of the dynamical equations.However, depending on δ it does not characterizes thedynamics of the system, since it becomes an unstablefixed point. In order to highligh it we show in Fig.(11)-(upper panels) the Jacobian spectrum for such steadystates. While it is an attractor for δ < /
2, it becomesa repulsor for δ > / δ = 1 / • Phase II: for couplings 1 / < δ < ω/κ < / m = ω/κm = ( − (cid:112) − ω/κ ) ) / m = − (cid:114) C − ω/κ ) −
43 ( m ) m = m m = − m / √ m i = 0 , i = 1 , , . (53) Figure 11: Jacobian spectrum for the steady state of phasesI and II along the line with ω/κ = 1 / δ . In theupper panels we show the real and imaginary terms of theJacobian eigenvalues for the steady states of phase I - Eq.(52)In the bottom panel we show the spectrum for phase II -Eq.(53). Except for δ = 1 /
2, in both cases the spectrum hasa four-fold degeneracy of eigenvalues with zero real part.
We show in Fig.(10) the occupation numbers along thisphase. The steady states do not depend on the coupling δ and the occupation numbers for energy levels 1 and 3are now equal, with a nonzero occupation for energy level2. Similar to phase I, the conditions of Eq.(53) is valid asa steady state for a larger range in the phase diagram,specifically, it is a solution of the dynamical equations ∀ δ with ω/κ < /
3, though not always stable. In Fig.(11)-(bottom panels) we show the Jacobian spectrum for suchsteady states. We see that it corresponds to a repulsor(unstable steady state) for δ < /
2, while becoming anattractor for δ > /
2. At the transition point δ = 1 / • Phase III: for couplings 1 / < δ < ω/κ (cid:39) / ω/κ > SU (2)case, but also in a region of the strong dissipative regime(2 / (cid:47) ω/κ < • Critical line: for coupling δ = 1 / ∀ ω/κ we ob-serve a very peculiar behavior. As discussed, the steadystates of Eqs.(52) are solutions of the dynamical equa-tions also for δ = 1 /
2. In this case, however, its Jaco-bian spectrum has only imaginary terms, correspondingto center steady states. It indicates that for a given initialstate close to these steady states, the dynamics shouldcorrespond to closed orbits, typical of the BTC’s foundin the collective d = 2 system. We find that this is in-deed the case. We show in Fig.(12)-(upper left panel) the4 t n i ( t ) n ( t ) n ( t ) n ( t ) t -0.2-0.100.10.2 Λ L y a p j t n i ( t ) n ( t ) n ( t ) n ( t ) t -0.2-0.100.10.2 Λ L y a p j m a x j Λ L y a p j Figure 12: (Left panels)
We show the dynamics of theoccuppation number for different initial conditions, with sys-tem parameters δ = 1 / ω/κ = 0 .
2: (upper left) for initialstates with n (0) = 1 (continuous curves) and n (0) = 1 (dot-ted curves); (bottom left) for initial states given by Eq.(52)with δ = 0 .
45 (continuous curves) and δ = 0 .
49 (dottedcurves). (Right panels)
Lyapunov spectrum obtained fromBenettin’s approach for the same system parameters: (upperright) using initial state n (0) = 1; (bottom right) for initialstate given by Eq.(52) with δ = 0 .
45. We show in the insetthe largest Lyapunov exponent decaying for long times, in alog-log scale. In both cases we used a time step dt = 10 − forour numerical simulation. . dynamics for two initial states close to the steady statesof Eq.(52) for δ = 1 /
2. Both show closed, but different,orbits.Considering the case of an initial condition far fromthese steady states, thus out of a Jacobian linear stabil-ity approach, we see now the existence of limit-cycles.We show in Fig.(12)-(bottom left panel) the dynamicsfor initial states with n i (0) = 1, for i = 3 or 2, bothfeaturing the same limit-cycle dynamics. We also ob-served that other different initial states may lead to evenfurther different limit-cycles. We remark that for largercoherent coupling ω/κ (cid:39) / X. CONCLUSION
In this work we studied boundary time crystals incollective d = 2 , d = 2 presented in [1] and extendedthe analysis of its phase diagram, obtaining the full setof steady states combining analytical and algebraically(quasi-analytically) approaches, and further studying itsJacobian stability. The existence of BTC in the model isseen to be directly related to the presence of center fixedpoints. Moreover, we obtained analytically the effectsof a specific Z symmetry breaking Hamiltonian term tothe dynamics, showing that BTC’s are destroyed by sucha perturbation in the model (see also [28] for a similarresult in a p , q interacting spin model).The analysis of the collective d = 4-level system, com-posed of a pair of collective interacting 2-level systems,showed even more fruitful. The model supports more ro-bust forms of BTC’s, from limit-cycles to period doublingbifurcations leading to chaos. The BTC is robust to Z symmetry breaking Hamiltonian term in this case. Weobtained the orbit diagram of the model from its collec-tive magnetization and extracted its bifurcation ratio b n for a finite n . The bifurcation ratio for finite n was founddifferent from the Feigenbaum constant of the seminal lo-gistic map. A careful analysis for the ratio in the limit n → ∞ remains as an interesting perspective for a futurework.In the collective d = 3-level system we observed thatdepending on the competition between the two Lindbal-dians L and L the model supports static steady statescharacterized by a “dark level”, limit-cycle dynamics ora peculiar dynamical phase at the critical line δ = 1 / d -level systemscan support different forms of BTC’s with richer proper-ties. It would be interesting e.g. to explore larger d ’s andits implications to these phases, possible applications aswell as a throughout analysis on the role of the globaland dynamical symmetries for such models. Acknowledgements
We acknowledge enlightening discussions with R.Fazio. F.I. acknowledges the financial support of theBrazilian funding agencies National Council for Sci-entific and Technological Development—CNPq (Grant5No. 308205 / / . / [1] F. Iemini, A. Russomanno, J. Keeling, M. Schir`o,M. Dalmonte, and R. Fazio, Phys. Rev. Lett. ,035301 (2018), URL https://link.aps.org/doi/10.1103/PhysRevLett.121.035301 .[2] N. Goldenfeld (1992).[3] S. Sachdev, Cambridge University Press (2000).[4] F. Wilczek, Phys. Rev. Lett. , 160401 (2012), URL https://link.aps.org/doi/10.1103/PhysRevLett.109.160401 .[5] T. Li, Z.-X. Gong, Z.-Q. Yin, H. T. Quan, X. Yin,P. Zhang, L.-M. Duan, and X. Zhang, Phys. Rev. Lett. , 163001 (2012), URL https://link.aps.org/doi/10.1103/PhysRevLett.109.163001 .[6] P. Bruno, Phys. Rev. Lett. , 118901 (2013), URL https://link.aps.org/doi/10.1103/PhysRevLett.110.118901 .[7] P. Bruno, Phys. Rev. Lett. , 029301 (2013), URL https://link.aps.org/doi/10.1103/PhysRevLett.111.029301 .[8] P. Bruno, Phys. Rev. Lett. , 070402 (2013), URL https://link.aps.org/doi/10.1103/PhysRevLett.111.070402 .[9] Nozi`eres, Philippe, EPL , 57008 (2013), URL https://doi.org/10.1209/0295-5075/103/57008 .[10] G. Volovik, Jetp Lett. , 491–495 (2013), URL https://doi.org/10.1134/S0021364013210133 .[11] H. Watanabe and M. Oshikawa, Phys. Rev. Lett. ,251603 (2015), URL https://link.aps.org/doi/10.1103/PhysRevLett.114.251603 .[12] A. Russomanno, F. Iemini, M. Dalmonte, and R. Fazio,Phys. Rev. B , 214307 (2017), URL https://link.aps.org/doi/10.1103/PhysRevB.95.214307 .[13] F. M. Surace, A. Russomanno, M. Dalmonte, A. Silva,R. Fazio, and F. Iemini, Phys. Rev. B , 104303 (2019),URL https://link.aps.org/doi/10.1103/PhysRevB.99.104303 .[14] K. Sacha, Phys. Rev. A , 033617 (2015), URL https://link.aps.org/doi/10.1103/PhysRevA.91.033617 .[15] A. Syrwid, J. Zakrzewski, and K. Sacha, Phys. Rev. Lett. , 250602 (2017), URL https://link.aps.org/doi/10.1103/PhysRevLett.119.250602 .[16] N. V. Prokof’ev and B. V. Svistunov, Journal of Exper-imental and Theoretical Physics , 860 (2018), URL https://doi.org/10.1134/S1063776118110092 .[17] D. V. Else, B. Bauer, and C. Nayak, Phys. Rev. Lett. , 090402 (2016), URL https://link.aps.org/doi/10.1103/PhysRevLett.117.090402 .[18] V. Khemani, A. Lazarides, R. Moessner, and S. L.Sondhi, Phys. Rev. Lett. , 250401 (2016), URL https://link.aps.org/doi/10.1103/PhysRevLett.116.250401 .[19] C. W. von Keyserlingk, V. Khemani, and S. L. Sondhi,Phys. Rev. B , 085112 (2016), URL https://link.aps.org/doi/10.1103/PhysRevB.94.085112 .[20] N. Y. Yao, A. C. Potter, I.-D. Potirniche, and A. Vish-wanath, Phys. Rev. Lett. , 030401 (2017), URL https://link.aps.org/doi/10.1103/PhysRevLett.118.030401 . [21] W. W. Ho, S. Choi, M. D. Lukin, and D. A. Abanin,Phys. Rev. Lett. , 010602 (2017), URL https://link.aps.org/doi/10.1103/PhysRevLett.119.010602 .[22] B. Huang, Y.-H. Wu, and W. V. Liu, Phys. Rev. Lett. , 110603 (2018), URL https://link.aps.org/doi/10.1103/PhysRevLett.120.110603 .[23] A. Russomanno, S. Notarnicola, F. M. Surace, R. Fazio,M. Dalmonte, and M. Heyl, Phys. Rev. Research ,012003 (2020), URL https://link.aps.org/doi/10.1103/PhysRevResearch.2.012003 .[24] R. Khasseh, R. Fazio, S. Ruffo, and A. Russomanno,Phys. Rev. Lett. , 184301 (2019), URL https://link.aps.org/doi/10.1103/PhysRevLett.123.184301 .[25] N. Y. Yao, C. Nayak, L. Balents, and M. P. Zaletel, Na-ture Physics , 438 (2020), URL https://doi.org/10.1038/s41567-019-0782-3 .[26] R. Hurtado-Guti´errez, C. F. Carollo, P´erez-Espigares,and P. Hurtado (2020), URL https://arxiv.org/abs/1912.02733 .[27] P. Wang and R. Fazio, Phys. Rev. A , 013306 (2021),URL https://link.aps.org/doi/10.1103/PhysRevA.103.013306 .[28] G. Piccitto, M. Wauters, F. Nori, and N. Shammah(2021), URL arXiv:2101.05710 .[29] A. Riera-Campeny, M. Moreno-Cardoner, and A. San-pera, Quantum , 270 (2020), ISSN 2521-327X, URL https://doi.org/10.22331/q-2020-05-25-270 .[30] A. Lazarides, S. Roy, F. Piazza, and R. Moessner, Phys.Rev. Research , 022002 (2020), URL https://link.aps.org/doi/10.1103/PhysRevResearch.2.022002 .[31] C. Lled´o and M. H. Szyma´nska, New Journal of Physics , 075002 (2020), URL https://doi.org/10.1088/1367-2630/ab9ae3 .[32] K. Seibold, R. Rota, and V. Savona, Phys. Rev. A , 033839 (2020), URL https://link.aps.org/doi/10.1103/PhysRevA.101.033839 .[33] J. e. a. Zhang, Nature (2017), URL https://doi.org/10.1038/nature21413 .[34] S. e. a. Choi, Nature (2017), URL https://doi.org/10.1038/nature21426 .[35] J. Rovny, R. L. Blum, and S. E. Barrett, Phys. Rev. Lett. , 180603 (2018), URL https://link.aps.org/doi/10.1103/PhysRevLett.120.180603 .[36] J. Rovny, R. L. Blum, and S. E. Barrett, Phys. Rev.B , 184301 (2018), URL https://link.aps.org/doi/10.1103/PhysRevB.97.184301 .[37] S. Pal, N. Nishad, T. S. Mahesh, and G. J. Sreejith, Phys.Rev. Lett. , 180602 (2018), URL https://link.aps.org/doi/10.1103/PhysRevLett.120.180602 .[38] J. Smits, L. Liao, H. T. C. Stoof, and P. van der Straten,Phys. Rev. Lett. , 185301 (2018), URL https://link.aps.org/doi/10.1103/PhysRevLett.121.185301 .[39] K. Sacha and J. Zakrzewski, Reports on Progress inPhysics , 016401 (2017), URL https://doi.org/10.1088%2F1361-6633%2Faa8b38 .[40] D. V. Else, C. Monroe, C. Nayak, and N. Y.Yao, Annual Review of Condensed Matter Physics , 467 (2020), URL https://doi.org/10.1146/ annurev-conmatphys-031119-050658 .[41] M. P. Estarellas, T. Osada, V. M. Basti-das, B. Renoust, K. Sanaka, W. J. Munro,and K. Nemoto, Science Advances (2020),https://advances.sciencemag.org/content/6/42/eaay8892.full.pdf,URL https://advances.sciencemag.org/content/6/42/eaay8892 .[42] F. Minganti, I. I. Arkhipov, A. Miranowicz, and F. Nori(2020), URL arXiv:2008.08075 .[43] X. Yang and Z. Cai, Phys. Rev. Lett. ,020602 (2021), URL https://link.aps.org/doi/10.1103/PhysRevLett.126.020602 .[44] S. Autti, P. J. Heikkinen, J. T. Makinen, G. E. Volovik,V. V. Zavjalov, and V. B. Eltsov, Nature Materials , 171–174 (2020), URL https://doi.org/10.1038/s41563-020-0780-y .[45] R. W. Bomantara and J. Gong, Phys. Rev. Lett. ,230405 (2018), URL https://link.aps.org/doi/10.1103/PhysRevLett.120.230405 .[46] M. Nielsen and I. Chuang,Quantum Computation and Quantum Information,Cambridge Series on Information and the NaturalSciences (Cambridge University Press, 2000), ISBN9780521635035, URL https://books.google.com.br/books?id=65FqEKQOfP8C .[47] J. Hannukainen and J. Larson, Phys. Rev. A ,042113 (2018), URL https://link.aps.org/doi/10.1103/PhysRevA.98.042113 .[48] D. F. Walls, P. D. Drummond, S. S. Hassan,and H. J. Carmichael, Progress of Theoreti-cal Physics Supplement , 307 (1978), ISSN 0375-9687, https://academic.oup.com/ptps/article-pdf/doi/10.1143/PTPS.64.307/5292058/64-307.pdf,URL https://doi.org/10.1143/PTPS.64.307 .[49] P. Drummond and H. Carmichael, Optics Commu-nications , 160 (1978), ISSN 0030-4018, URL .[50] R. Puri and S. Lawande, Physics Letters A , 200 (1979), ISSN 0375-9601, URL .[51] D. F. Walls, Journal of Physics B: Atomic and Molecu-lar Physics , 2001 (1980), URL https://doi.org/10.1088/0022-3700/13/10/008 .[52] A. Pikovsky and A. Politi, Cambridge Univer-sity Press (2015), URL http://dx.doi.org/10.1017/CBO9781139343473 .[53] A. C. Y. Li, P. F., and K. Jens, Scientific Reports , 4887(2014), URL https://doi.org/10.1038/srep04887 .[54] (A detailed analysis of the initial state dependence is aninteresting topic. It is, however, beyond the scope of thiswork).[55] G. Benettin, L. Galgani, and J.-M. Strelcyn, Phys. Rev.A , 2338 (1976), URL https://link.aps.org/doi/10.1103/PhysRevA.14.2338 .[56] G. Benettin, L. Galgani, A. Giorgilli, and J.-M. Strel-cyn, Meccanica , 21 (1980), URL https://doi.org/10.1007/BF02128237 .[57] M. Sandri, The Mathematica Journal6