Quantum Neimark-Sacker bifurcation

Recently, it has been demonstrated that asymptotic states of open quantum system can undergo qualitative changes resembling pitchfork, saddle-node, and period doubling classical bifurcations. Here, making use of the periodically modulated open quantum dimer model, we report and investigate a quantum Neimark-Sacker bifurcation. Its classical counterpart is the birth of a torus (an invariant curve in the Poincaré section) due to instability of a limit cycle (fixed point of the Poincaré map). The quantum system exhibits a transition from unimodal to bagel shaped stroboscopic distributions, as for Husimi representation, as for observables. The spectral properties of Floquet map experience changes reminiscent of the classical case, a pair of complex conjugated eigenvalues approaching a unit circle. Quantum Monte-Carlo wave function unraveling of the Lindblad master equation yields dynamics of single trajectories on “quantumtorus” and allows for quantifying it by rotation number. The bifurcation is sensitive to the number of quantum particles that can also be regarded as a control parameter.

bagel-shaped form -a section of "quantum torus" -for the boson interaction strength close to the bifurcation value for the mean-field model. Importantly, the same transformation is observed in the stroboscopic distribution of observables, obtained by the Monte Carlo wave-function stochastic unraveling of the Lindblad equation, the method, especially relevant in the context of quantum optics and cavity systems [38][39][40][41] . Dynamics of such individual quantum trajectories on "quntum torus" is suitably characterized by rotation number. Similar to the classical case, rotation numbers close to rational correspond to almost "periodic" multi-modal stroboscopic distributions. Finally, we demonstrate that the bifurcation is system size dependent, which plays a role of another bifurcation parameter.
The paper is organized as follows. In Section 1 we describe the quantum model, its nonlinear mean-field approximation, and numerical methods. Section 2 contains the main results. Section 3 gives a summary and outlook.

Model and Methods
Within the Markovian approximation framework (which assumes weak coupling to environment), the evolution of an open quantum system can be described by the Lindblad master equation 21,23 , where the first term in the r.h.s. captures the unitary evolution, and the second term describes the action of environment. We consider a system of N indistinguishable interacting bosons, that hop between the sites of a periodically rocked dimer. This model is a popular theoretical testbed [42][43][44] , recently implemented in experiments [45][46][47] , known to exhibit regular and chaotic regimes [27][28][29]33 . Its unitary dynamics is governed by the Hamiltonian Here, J denotes the tunneling amplitude, U is the interaction strength, and ε t ( ) presents a periodical modulation of the on-site potentials. It is chosen as π Ω = T 2 / , so that the amplitude A is the dynamic energy offset between the two sites. Here b j and † b j are the annihilation and creation operators of a particle at site j, while = † n b b j j j . The system Hilbert space has dimension + N 1 and can be spanned with + N 1 Fock basis vectors, labeled by the number of bosons on the first site, | 〉 n { }, = … n N 0, , . Thus, the size of the model is controlled by the total number of bosons.
The dissipative term involves a single experimentally relevant jump operator [48][49][50] : which attempts to 'synchronize' the dynamics on the two sites by constantly recycling anti-symmetric out-phase modes into symmetric in-phase ones. The dissipative coupling constant γ is taken to be time-independent. Throughout the paper we set = J 1, γ = . 0 1, = .
A 3 4, π = T 2 , and vary U and N. The computational analysis makes use of two methods to evolve the system numerically. First, we implement propagation of Eq. (1) by the fourth-order Runge-Kutta scheme. The evolution converges to a unique asymptotic solution, which in case of periodic modulation is a stable periodic trajectory (or a fixed point of the corresponding stroboscopic map, ) F , = … m 0, 1, 2, ). For the particular system Eqs. (2) and (4), we set the integration time step of ⋅ − T 5 10 4 , and leave at least 100 modulation periods for transients to fade out. The density matrix ρ of the system with N bosons can be visualized on the Bloch sphere by plotting the Husimi distribution ϑ φ p( , ), which can be obtained by projecting ρ on the set of the generalized SU(2) coherent states 51,52 The corresponding nonlinear mean-field equations for the θ ϕ ( , ) phase variables are obtained by writing the master equation in terms of the spin operators = www.nature.com/scientificreports www.nature.com/scientificreports/ where proportional to γ terms of lower order in N are neglected. The quantity = + + S S S S x y z 2 2 2 2 is a constant of motion, so the mean-field evolution is restricted to the surface of a Bloch sphere, φ ϑ = S S S ( , , ) [ cos( )sin( ), The corresponding particle number in the first site is then recovered as . This nonlinear dynamical system plays a reference role in the bifurcation analysis further. In numerical experiments we record the consecutive stroboscopic values, 0, 1, 2, , and construct the histograms normalized to the largest value 1 for each parameter set, complemented by Poincaré sections for the stroboscopic map.
Second, we employ the Monte-Carlo wave function (also "quantum trajectory" or "quantum jump") method 38,39 to unravel Lindblad master Eq. (1) into an ensemble of quantum trajectories. It recasts the evolution of the model system into the ensemble of systems described by wave functions, ψ t ( ) r , = … r M 1, 2, , r , governed by an effective non-Hermitian Hamiltonian, H . This Hamiltonian incorporates the dissipative operator V, which is responsible for the decay of the norm, When the norm drops below a threshold, that is randomly chosen each time, a quantum jump is simulated, such that the wave function is transformed according to ψ ψ → V and then normalized 20 . The density matrix can then be sampled from a set of M r realizations as ρ time for relaxation towards an asymptotic state, and following the dynamics for up to = t T 10 3 . Thus we produce stroboscopic plots for the expectation values of the two observables, the number of particles on the left site of the dimer, n(t), and the energy, e(t), Quantum trajectories allow for an insight to dynamics in the asymptotic regime -on quantum attractorbeyond the stationary Husimi projection picture.

Results
We first investigate the nonlinear mean-field equations, Eq. (7), treating interaction strength U as a bifurcation parameter. There one observes the emergence of an invariant curve for a stroboscopic map (  1) )} by the Neimark-Sacker bifurcation at ∼ . U 0 11 and its succession by period-6 cycle at ∼ . U 0 18 (Figs. 1(a) and 2). Further it develops into a chaotic attractor, that ultimately disappears through a crisis, when the stable fixed point is recovered.
We interrogate whether quantum equations manifest an analogue of the classical Neimark-Sacker bifurcation and study its properties.
A sufficiently large number of particles, = N 500, should bring the system close to the classical limit. The one-parameter bifurcation diagram displays the probabilities to find a given number of particles on the left site of a dimer, taken at stroboscopic times mT, = … m 0, 1, 2, , numerically given by the diagonal elements of the density matrix, ρ mT ( ) n n , . The qualitative structure of the bifurcation diagram obtained from the mean-field model is reproduced in the quantum case quite well, as for the birth of a torus, as for the onset of chaos and ultimate recovery of a stable fixed point, cf. Fig. 1(b). However, even for = N 500 some details differ, like the quantum bifurcation within ∈ .
. U [0 6, 0 7], that does not have a counterpart in the mean-field model. The discrepancy stems from truncation of double and higher order correlators, 1 , in the mean-field approximation, retaining only expectation values, cf. Eqs. (6 and 7). We conjecture that by going to a higher-order mean-field approximation, one would obtain a higher-dimensional nonlinear system, where the corresponding classical bifurcation would emerge.
Next, we demonstrate the correspondence between the asymptotic states at stroboscopic Poincaré sections (for mean-field model) on the plane ϑ φ { , } and Husimi distributions for the quantum system, Fig. 2. Indeed, the projection of the quantum attractor onto the classical phase space reveals the emergence of a bagel-shaped distribution after the invariant curve in the Poincaré section of the mean-field equations. Moreover, one observes the formation of a multi-period orbit on the quantum torus later on, cf. Fig. 2(d), the scenario typical of classical systems.
There are, however, quantitative differences: the bagel is already present at = . U 0 1 for the quantum model with = N 500, cf. Fig. 2(a), while the mean-field model still has a fixed point; the size of the bagel is slightly greater than the invariant curve in the classical case, cf. Fig. 2(b,c); the formation of a periodic orbit in the quantum case occurs at lower U, cf. Fig. 2(d). (2019) 9:17932 | https://doi.org/10.1038/s41598-019-53526-2 www.nature.com/scientificreports www.nature.com/scientificreports/ We further investigated, whether quantum Neimark-Sacker bifurcation can be identified directly in quantum observables. Sampling attractor with quantum trajectories (cf. Section 1), we reconstruct histograms for stroboscopic observables n mT e mT ( ), ( ), Eq. (10), and again witness the emergence of bagel-shaped distributions following quantum bifurcation, cf. Fig. 3 and compare to Fig. 2.
Quantum trajectory unraveling allows to get a deeper insight into the dynamics on quantum torus. That is, one can define the phase θ m for each pair of stroboscopic observables, n mT e mT ( ), ( ), as a polar angle, the origin placed at the center of mass of the stroboscopic 2D histogram, which is normalized in both dimensions to the same interval − [ 1, 1] (see Fig. 4, inset), and thus calculate the instantaneous rotation (or winding) number  In classical dynamics, the time-averaged ω discriminates two cases. For rational ω = p q / , ∈ p q ,  , the invariant torus contains a stable period-q orbit that is observed as an asymptotic solution, while for irrational ω the trajectories cover the torus densely. The regimes interchange with system parameters. 0 58, different from the period doubling value 1/2, the bifurcation present for the other parameter values 28,33 . Our case is analogous to irrational rotation number for classical torus, so that quantum trajectories densely cover the bagel. The average rotation number increases with interaction strength, and becomes rational ω = 3/5 at ≈ . U 0 15. Then the stroboscopic distributions obtain a clear footprint of the corresponding period-5 structure (Figs. 2(c,d) and 3(c,d)), as in the classical case.
The bifurcation also affects the spectral properties of the system. This can be demonstrated by calculating and diagonalizing the Floquet map where T T is time-ordering operator, that describes evolution of the density operator over a period of modulation under Eq. (1) 27 . The largest eigenvalue of its spectrum, μ { } k , = … + k N 1, , ( 1) 2 , is always unity, μ = 1 1 , the rest are inside the unit circle, as per the dissipative nature of L. We follow the next to largest eigenvalues and find a conjugated pair that approaches the unit circle at the point of quantum bifurcation, μ ≈ θ ± e i 2,3 0 (Fig. 5, inset). Its complex phases are consistent with the rotation number, θ πω ≈ 2 0 . Repeating the treatment of ref. 33 , one obtains that the two-time correlation of an observable in the asymptotic regime is incommensurate with the time-scale of modulation. Note that the spectral gap μ − | | 1 2,3 decreases with N, yielding the classical conjugate multipliers of the limit cycle at the Neimark-Sacker bifurcation point 2 in the infinite-size limit, → ∞ N . It also gives an estimate of relaxation time from a random initial condition to the asymptotic state, , which can become very large with N, it follows. Finally, we investigate the dependence on the number of bosons, N. The classical limit is formally obtained for → ∞ N , but the traits of classical Neimark-Sacker bifurcation are reproduced deep in the quantum regime, as soon as ∼ N 25. In fact, the number of particles, N, can be considered as a bifurcation parameter itself. For instance, one can explore the case = .
U 0 1125, when an invariant curve is already present in the mean-field model, cf. Figs. 2(b) and 3(b). However, the Husimi distribution for the quantum dimer with the relatively small number of particles, = N 50, is still unimodal, cf. Fig. 6(a). Increasing N one observes a transformation to a bagel shape ( Fig. 6(b-d)), i.e. Neimark-Sacker bifurcation with the system size N.
In classical nonlinear systems the size of the invariant curve generically scales as a square root above the bifurcation point 2 . For the quantum bifurcation we define the diameter of the bagel D in the Husimi distribution as the distance between two maxima in the section ϕ π = /2. When the distribution is unimodal, we set = D 0. The resulting curves D(U) confirm a pronounced dependence on the number of particles, in particular, for the bifurcation point, cf.  U 0 12, = N 50; red points mark μ 2,3 .

conclusions
We found and investigated the quantum counterpart of the classical Neimark-Sacker bifurcation in an open periodically modulated quantum dimer. The classical bifurcation in a dissipative nonlinear system is a birth of torus from a limit cycle (or an invariant curve from a fixed point for a time-discrete map). Concurrently, a conjugate pair of complex Floquet multipliers of a periodic trajectory crosses the unit circle. The quantum Neimark-Sacker bifurcation is manifested by emergence of bagel-shaped stroboscopic distributions, as for Husimi projection on the classical phase space, as for quantum observables. The bifurcation is also seen in the spectral properties of the corresponding Floquet map, as a pair of its conjugate eigenvalues approaching the unit circle. Quantum trajectories technique allows for unraveling the dynamics on the "quantum torus", and a suitably generalized rotation number on it. In analogy to the classical case, irrational and rational rotation numbers yield bagel (dense torus) and multi-modal (periodic orbit on torus) distributions. The quantum bifurcation depends on the system size -the number of bosons, and more sensitively for fewer particles, away from the mean-field limit. As such the number of particles is a bifurcation parameter.
We expect that the bifurcation can be found in the other open non-equilibrium quantum setups, to name a modulated open Dicke model, where its classical counterpart has been reported 54 . Thereby, cavity and circuit QED systems appear the first candidates for its experimental observation. U 0 1125.