The driven-dissipative Bose-Hubbard dimer: phase diagram and chaos

We present the phase diagram of the mean-field driven-dissipative Bose-Hubbard dimer model. For a dimer with repulsive on-site interactions ($U>0$) and coherent driving we prove that $\mathbb{Z}_2$-symmetry breaking, via pitchfork bifurcations with sizable extensions of the asymmetric solutions, require a negative tunneling parameter ($J<0$). In addition, we show that the model exhibits deterministic dissipative chaos. The chaotic attractor emerges from a Shilnikov mechanism of a periodic orbit born in a Hopf bifurcation and, depending on its symmetry properties, it is either localized or not.

The Bose-Hubbard model is a celebrated fundamental quantum mechanical model that accounts for boson dynamics in a lattice [1]. It successfully describes the interplay between the hopping of particles between neighboring sites of the lattice (with rate J) and on-site interactions. Such interactions appear as multi-boson terms in the Hamiltonian with interaction energy U . Importantly, this model accurately explains the superfluid to Mott insulator phase transitions, that has been experimentally demonstrated in ultracold atomic lattices [2]. A minimal building block in this context is the so-called Bose Hubbard dimer, consisting of only two interacting sites, also known as the bosonic Josephson junction [3,4]. Futhermore, the Bose-Hubbard dimer lies at the basis of a number of striking phenomena such as the Josephson effect, self-trapping and symmetry breaking [5][6][7].
In recent years there has been a growing interest in understanding open quantum systems, where the bosons can be added and destroyed by means of external driving and dissipation mechanisms [8][9][10][11]. In this context, photonic systems have attracted much attention since photons in optical cavities can be injected through an external driving laser, and dissipation comes in as a natural consequence of optical cavity losses. The drivendissipative Bose-Hubbard dimer has been realized in a number of experimental systems, such as semiconductor microcavities [12,13], superconducting circuits [14] and photonic crystals [15], where the interactions take the form of Kerr-type optical non-linearities. In many regards, these optical systems constitute outstanding platforms for studying many-body phenomena in open quantum systems [16]. Among them, dissipative phase transitions [17,18] are an especially exciting open topic, because they provide a conceptual basis for the understanding and the prediction of new collective states, both steady and dynamical ones, with the latter accounting for collective coherent oscillations.
As is well known, phase transitions are characterized by critical phenomena, which can emerge in the thermodynamic limit, i.e., with a large photon number in optical cavities. Remarkably, even single-mode cavities -i.e., with no spatial degrees of freedom -may display phase transitions including optical bistability (which is of first order) [18,19] and the emergence of an oscillation threshold in, e.g., two-photon pumped Kerr resonators [20] or laser devices [21] (which is a second-order phase transition). As pointed out in Ref. [20], driven-dissipative phase transitions are strongly related to mean-field semiclassical solutions in such a thermodynamic limit: these are known as bifurcations of a classical vector field from a non-linear dynamical point of view; see, e.g., [22,23].
Recent examples of complex dynamics found in phase diagrams of driven-dissipative nonlinear cavities include oscillating phases [24] and exotic attractors [25]. Moreover, Lorenz-type chaos has already been found in the Gross-Pitaevsky equation accounting for incoherent pump of a double potential well Bose-Einstein condensate [26]. However, to our knowledge, dissipative chaos has not been predicted so far in the coherently driven regime, which is the one that is well described by the driven-dissipative Bose-Hubbard dimer. In this work, we show that dissipative chaos exist in this model and that it is intimately related to Shilnikov homoclinic bifurcations. In addition, we present the phase diagram as a result of the comprehensive analysis of local bifurcations of steady states and their symmetry properties.
The mean-field approximation of the Bose-Hubbard model is (see, e.g., [27]) where α 1 , α 2 ∈ C are the electric field envelopes in each optical cavity, which are linearly coupled by the tunneling parameter J. Both cavities are coherently driven by a field with amplitude F and detuning ∆ = ω p − ω c , where ω p and ω c are the driving and cavity frequencies, respectively. The cavities have a Kerr-type nonlinearity characterised by U , which is assumed to be fast compared to the dissipation rate γ. We apply to system (1) the coordinate transformation i is the complex conjugate of α i . With time rescaled to τ = γt/2, we obtain the non-dimensionalized Bose-Hubbard model as the vector field for the (rescaled) envelopes A, B ∈ C. Here are the (rescaled) coupling, drive and detuning parameters, respectively.
Since system (2) is invariant under the transformation all results for positive U directly translate to those for negative U ; hence, we consider U > 0 and sign(U ) = 1 in (2) from now. Note that, although U is fixed by the material of the cavity, both the magnitude and sign of κ can be changed by means of photonic design in some particular geometries, for example, with potential barrier engineering for the case of photonic crystal nanocavities [15,28]. Because δ can be altered at will during the experiment by changing the frequency of the driving laser, states of negative interaction energy in an otherwise positive-U system can be assessed via the parameter transformations κ → −κ and δ → −δ. In addition, system (2) is invariant under the phasespace mirror symmetry of exchanging the two cavities, that is, of swapping A and B, which gives rise to Z 2equivariance of solutions [23]. Hence, solutions (equilibria, periodic orbits, trajectories) can be split into two major groups: symmetric and asymmetric ones. For a symmetric solution the field intensities and phases in both cavities are the same. The set of all symmetric solutions forms the symmetry subspace given by A = B, which is an invariant subspace of system (2). An asymmetric solution, on the other hand, is one where the intensity and/or the phase in both cavities differ, that is, A = B. Note that asymmetric solutions often come in pairs, one being the mirror image of the other under swapping A and B.
We start our analysis of system (2) by characterizing the equilibria and their bifurcations as the drive parameter f increases. Figure 1 shows phase diagrams of |A| 2 | A| 2 6.0 against f for three different values of δ for the two cases of a negative and positive coupling parameter κ in the left and right columns, respectively. The corresponding phase diagrams of |B| 2 against f are the same due to Z 2 -equivariance. All bifurcations and branches of solutions have been computed with the numerical continuation package Auto07p [29]. When f = 0 there is the stable equilibrium given by A = B = 0, which gives rise to a monotone branch of stable symmetric equilibria for any f provided the detuning δ is sufficiently large. However, when δ = 0 and for increasingly negative δ, one finds bifurcations on the branch of symmetric equilibria.
For negative κ, as in Fig. 1a-c for κ = −3.5, we observe for δ = 0 an interval of bistability that is delimited by two saddle-node bifurcations S of symmetric equilibria [ Fig. 1a]. In this interval there are three symmetric equilibria that form a hysteresis loop: a stable equilibrium corresponding to low intensity in both cavities, an intermediate saddle equilibrium, and a stable equilibrium with higher intensity in both cavities. As δ decreases, the interval of bistability increases in size. While bistability can be encountered in one-cavity systems, spontaneous Z 2 -symmetry breaking requires coupled cavities with mirror symmetry; this is the case for the Bose Hubbard model (2) and spontaneous symmetry breaking phase transitions are known to exist in the form of pitchfork bifurcations [27,30]. As Fig. 1b for δ = −6 shows, we find two such symmetry breakings at the points labelled P on the low intensity branch. The left pitchfork bifurcation point P is subcritical and gives rise to a pair of unstable asymmetric equilibria, while the right point P is supercritical and gives rise to a pair of stable asymmetric equilibria (one with |A| 2 < |B| 2 and the other with |A| 2 > |B| 2 ); these two branches come together at a pair of saddle-node bifurcations of asymmetric equilibria denoted S * . These additional bifurcations create a small f -interval, between S * and the left point P, with four stable equilibria: higher and lower intensity symmetric equilibria, and a pair of asymmetric (lowerintensity) equilibria [ Fig. 1b]. Notice that between the two points P there are still three stable objects: the two asymmetric equilibria and the higher-intensity symmetric equilibrium. For an even lower detuning, as in Fig. 1c for δ = −6.2, we find two pairs of supercritical Hopf bifurcations denoted H on the stable branch of asymmetric equilibria. They give rise to a pair of stable asymmetric periodic orbits, which are represented in the phase diagram by their maxima in |A| 2 . In the corresponding f -interval the asymmetric equilibria are now unstable. Note that the two asymmetric periodic orbits are new attractors that coexist with the stable higher intensity equilibrium.
The situation for κ > 0 is notably different: bistability of the symmetric state is, in fact, absent for all panels Fig. 1d-f. Instead, for δ = 0 we find symmetry breaking at two supercritical pitchfork bifurcations, which give rise to a pair of stable asymmetric equilibria in an quite large f -interval where the symmetric equilibrium is unstable [ Fig. 1d]. This makes this set of abnormal coupling conditions well suited for experimental implementations. Other driving configurations might also allow one to fulfill these conditions, but then the driving phase of the cavities needs to be adjusted properly such that the antisymmetric state can be linearly excited [27].
As δ is decreased, as in Fig. 1e for δ = −3.0, the pair of asymmetric stable equilibria that exists in this interval of symmetry breaking also exhibit saddle-node bifurcations S * that create bistability of asymmetric states and a pair of hysteresis loops. Perhaps surprisingly, there is an frange where the intensity in one of the cavities is nearly zero. As Fig. 1f for δ = −4.5 shows, a high negative detuning generates a pair of supercritical Hopf bifurcation points labelled H, which create an f -interval with pair of asymmetrical stable periodic solutions, while the asymmetric equilibria are unstable.
The bifurcations of equilibria identified in Fig. 1 can be followed in the additional parameter δ, yielding a phase diagram in the (f, δ)-parameter plane consisting of bifurcation curves that bound regions with qualitatively different behavior. Figure 2  alytical expressions obtained by solving for the equilibria and the respective bifurcation condition [31]. These calculations also show that the saddle-node bifurcation S in the symmetry subspace (brown line) only occurs if δ < − √ 3 − κ and 1.24081 < f . At the points (f, δ) ≈ 1.24081, − √ 3 − κ there is a cusp bifurcation CP on the curve S that corresponds to a sharp bound for bistability in the (f, δ)-plane [23]. The region of bistable symmetric equilibria bounded by S is shaded brown in both panels of Fig. 2; note the quite different position of this region for κ = −3.5 and for κ = 3.5. Additionally, we find that the curve of pitchfork bifurcation P (purple curve), which bounds the shaded region with asymmetric equilibria, has a maximum at δ = − √ 3 + κ; hence, symmetry breaking only occurs if δ < − √ 3 + κ, which is equivalent to the bound found in [32] and also consistent with [33]. As an unexpected consequence, the condition for symmetry breaking at P to be possible is that the detuning ω p − (ω c − κγ/2), with respect to the antisymmetric mode, must exceed √ 3/2 times the cavity linewidth. This is similar to the bistability condition, but the detuning is measured from the antisymmetric mode of the system, i.e., the one that cannot be linearly excited from the outside world. These bounds explain why the first phenomenon observed, as δ decreased, is bistability when κ < 0, while it is symmetry breaking when κ > 0; compare the two panels of Fig. 2 and see also Fig. 1. Notice that the curve P is tangent to the curve S at the saddle-node pitchfork point SP, to the right of the cusp point CP for κ < 0 and to the left of CP for κ > 0; moreover, for κ < 0 the curve P changes criticality at the point DP, from which the curve S * of saddle-node bifurcations of asymmetric equilibria emerges.
The curves SN * and H in Fig. 2 of saddle-node and Hopf bifurcations or asymmetric equilibria are computed numerically by continuation techniques. They bound regions of bistability and of stable asymmetric periodic solutions, respectively. Notice that for κ < 0 the curve H emanates from a saddle-node pitchfork point SP and is superctitical throughout, while for κ > 0 we find a change of criticality of H at a generalised Hopf bifurcation (labelled GH) [23]. The considerable amount of multistability of the Bose-Hubbard model (2) is represented by the overlap between the different shaded region of Fig. 2. For even lower negative values of δ than shown in Fig. 2, one can find up to nine equilibria, as pointed out in [30].
As δ is decreased for both signs of κ, the asymmetric stable periodic orbits created at the supercritical Hopf bifurcation H exhibit bifurcations scenarios to more complex and chaotic dynamics; the region of complex dynamics is indicated by grey shading in Fig. 2. Its exact boundary is formed by an intricate arrangement of accumulating bifurcation curves, including period-doubling and different types of homoclinic and heteroclinic bifurcations [31]. Notice that for κ < 0 periodic solutions become more complex almost immediately, while for κ > 0 we find a quite large region of stable periodic solutions. While a full discussion of the different scenarios is beyond the scope of this paper, we conclude by showing the existence of two types of Shilnikov-type chaos in the Bose-Hubbard model: asymmetric chaotic dynamics localized in one of the cavities, and symmetric chaotic dynamics with irregular switching between the two cavities.
At a homoclinic bifurcation a system of differential equations posesses a special trajectory, called a homoclinic orbit, that converges both in backward and forward time to a saddle equilibrium. For the special case that the equilibrium is a saddle focus (with a pair of complex conjugate eigenvalues), a celebrated result by Shilnikov [34] states that (under a certain eigenvalue conditions) infinitely many periodic orbits and chaotic behavior can be found nearby. This phenomenon is now referred to as a Shilnikov bifurcation. In the Z 2 -equivariant system (2) we find Shilnikov bifurcations both to asymmetric saddle foci and to symmetric saddle focus, as well as two types of nearby chaotic attractors. Figure 3 shows both cases of an asymmetric and of a symmetric Shilnikov bifurcation and associated chaotic attractors. For κ < 0 in Fig. 3a, there exist asymmetric Shilnikov homoclinic orbits to each of the pair of asym- metric equilibria, one either side of the symmetry line |A| 2 = |B| 2 . For nearby f as in Fig. 3b, we find a pair of chaotic asymmetric and localized attractors. The time traces in panels (c1,c2) show that this type of localized chaotic dynamics is characterized by a dominance in intensity of one of the two cavities for all time. For κ > 0, on the other hand, we find a pair of Shilnikov homoclinic orbits to a single symmetric saddle equilibrium [ Fig. 3d]. This has the consequence that the associated chaotic attractor in Fig. 3e is much larger and no longer localized as it crosses the line |A| 2 = |B| 2 . As the time traces in panels (f1,f2) clearly show, this type of chaotic dynamics is characterized by episodes of dominance in intensity of one of the two cavities for some time, with rapid and irregular switching events of the intensity to the other cavity. We remark that a local analysis of equilibria alone is not able to predict and explain this global switching behavior of the system. In conclusion, our systematic bifurcation analysis pro-duced a comprehensive phase diagram description of both equilibria and bifurcating non-equilibrium solutions of the driven-dissipative Bose-Hubbard dimer model. In particular, we were able to identify bifurcations to periodic solutions, and subsequent localized and nonlocalized chaotic behavior near homoclinic bifurcations of Shilnikov type. Crucial to the dynamics is the fact that the two cavities are coherently driven, which means that the system remains four-dimensional since there is no phase freedom as in active cavities. Interestingly, in a system with positive on-site interaction energy (U > 0), non-localized Shilnikov chaos can be observed for negative tunneling parameter (J < 0), in a bifurcation cascade whose details will be reported elsewhere. Our meanfield dynamical results thus pave the way for studying complex non-equilibrium orbits -including quantum dissipative chaos -in the quantum master equation of the open Bose-Hubbard dimer. This work has been partially funded by the "Investissements d'Avenir" program (Labex NanoSaclay, Grant No. ANR-10-LABX-0035) and the ANR UNIQ DS078.