Semiclassical bifurcations and quantum trajectories: a case study of the open Bose–Hubbard dimer

We consider the open two-site Bose–Hubbard dimer, a well-known quantum mechanical model that has been realised recently for photons in two coupled photonic crystal nanocavities. The system is described by a Lindblad master equation which, for large numbers of photons, gives rise to a limiting semiclassical model in the form of a four-dimensional vector field. From the situation where both sites trap the same amount of photons under symmetric pumping, one encounters a transition that involves symmetry breaking, the creation of periodic oscillations and multistability as the pump strength is increased. We show that the associated one-parameter bifurcation diagram of the semiclassical model captures the essence of statistical properties of computed quantum trajectories as the pump strength is increased. Even for small numbers of photons, the fingerprint of the semiclassical bifurcations can be recognised reliably in observables of quantum trajectories.


Introduction
Phase transitions describe a fundamental change in the behaviour of a system as its parameters are changed. Understanding why these systems exhibit different observable features in different physical scenarios and how they transition between them is a fundamental problem in many branches of physics. In quantum systems, the interest is in phase transitions, which occur near zero temperature and are driven by quantum fluctuations [20,26]. We are concerned here with an open quantum system, where particles of an ensemble may be lost to the environment. Open quantum systems constitute a realistic scenario from an experimental perspective, where lost particles can be measured to monitor the system as it evolves. The underlying non-unitary evolution can produce dissipation and decoherence not exhibited by their closed counterparts.
For both open and closed quantum systems, quantum fluctuations can become negligible in the thermodynamic limit of large particle numbers (e.g., atoms and/or photons). Their time evolution can be described by a so-called semiclassical model or mean-field approximation [42], which takes the form of a closed set of ordinary differential equations (ODEs) for associated averaged quantities. Such an approach has long been common in quantum optics to study systems with large numbers of photons, such as lasers, and it is the starting a e-mail: a.giraldo@auckland.ac.nz (corresponding author) point of our analysis here. The attractors of the semiclassical ODE, that is, its stable solutions, therefore, give insight into the observable behaviour of the underlying quantum system. In other words, phase transitions exhibited by the quantum system can be identified as changes in the stability of solutions of the limiting semiclassical model, which are examples of bifurcations. Thus, bifurcation analysis of the semiclassical ODEwith a combination of analytical and advanced numerical tools-allows one to systematically map out the attractors and phase transitions in the thermodynamical limit; see for example [21,32,33] as entry points to bifurcation theory and associated numerical methods. It has long been common in quantum optics to study systems with large numbers of photons, such as lasers, by means of their semiclassical descriptions.
We are interested here in the relationship between dynamics and bifurcations of the semiclassical ODE model and the observable behaviour of the quantum system when the number of particles of the considered system is relatively small. In this case, the system is far from the thermodynamic limit, and quantum fluctuations may be significant and cannot simply be neglected. Quantum simulations of the underlying Hamiltonian, which are only feasible computationally for small numbers of particles, can be used to investigate the predictive power of semiclassical models for novel types of quantum systems that operate with very few atoms and/or photons. It has been shown by means of quantum simulations that fingerprints of semiclassi- cal predictions can still be found in probabilistic features of quantum systems for surprisingly low numbers of particles [7,22,49,51]. As a recent example, the semiclassical limit of the unbalanced Dicke model [48] predicts parameter regimes of superradiant switching, quantum hysteresis, and oscillations, which have been observed in its quantum optical description [49].
This paper probes the connections between the semiclassical and quantum regimes for the specific example of the open two-site Bose-Hubbard dimer [14,27,31]. This well-known quantum mechanical model describes the dynamics of bosonic particles in a lattice, where the behaviour is determined by the interplay of the hopping rate of particles between lattice sites and their on-site interactions. Different experimental realisations of this system have been achieved in the form of semiconductor microcavities [1,34,46] and superconducting circuits [13,44]. Of particular interest to us is an optical realisation in the form of two lossy coupled photonic crystal resonators [16,24]. In this setup, good agreement between experimental measurements and the semiclassical description has been obtained in a regime up to moderate strength of the optical drive signal, in particular, concerning the observation of spontaneous symmetry breaking [18]. Recent theoretical results over a much wider range of parameters have uncovered a rich variety of behaviour of the semiclassical open Bose-Hubbard dimer, including chaotic behaviour with non-switching, regular switching and chaotic switching between the two sites; see [17] for more details. We focus here specifically on the comparison between the quantum and semiclassical descriptions of the two-site Bose-Hubbard model. Such comparisons have been performed before under the excitation of the anti-bonding mode [8], and under excitation of the bonding mode for the case of positive intermode coupling [5,37]. Here we consider the situation, where the bonding mode is only excited, while the intermode coupling of the cavities is negative. In this way, we extend earlier comparisons to the specific parameter constellations that most resembles the experimental setup in [16] and, to our knowledge, have not been studied before.
Our approach is as follows: we choose parameter values from [18] to obtain a bifurcation diagram of the semiclassical ODE model with different qualitative behaviour, specifically, symmetry breaking, periodic motion and multi-stability of asymmetric states. We then obtain temporal traces of the open Bose-Hubbard dimer by means of quantum trajectory simulations [6]. We then compare the two. By performing simulations of the quantum system for an increasing number of photons, we are able to observe the emergence of semiclassical attractors in a fully quantum realisation even with low photon numbers. Furthermore, we consider signatures for the existence of anti-bunching [39] and entanglement [40,43] in the quantum system as different bifurcations of the semiclassical approximation are encountered. Overall, our results showcase the role a semiclassical bifurcation analysis of a quantum system may play in providing a roadmap of interesting fundamental behaviour of interest, both from a theoretical and an experimental perspective. Specifically for the open two-site Bose-Hubbard dimer, this demonstrates that it might be feasible to investigate the quantum footprint of even quite complex dynamical behaviour, such as different types of chaotic switching [17].
The computations of equilibria and local bifurcations presented here are implemented in and performed with the pseudo-arclength continuation package Auto07p [10,11] and its extension HomCont [9]. The quantum trajectory simulations were carried out in the software package QuTiP [28,29]. Visualisation and postprocessing of the data are performed with Matlab R .

Quantum Hamiltonian and semiclassical ODE model
The Hamiltonian of the Bose-Hubbard model for two cavities takes the form [5]: whereâ † j is the creation operator for cavity j. The intermode coupling is represented by J, the frequency inside the cavities by ω c and the on-site energy by U . The last term of (1) represents a driving field with frequency ω p and amplitude F , which we consider equal for both cavities. In what follows, we work in a co-moving frame with respect of ω p ; in this way, the time dependency of the Hamiltonian (1) is dropped. Since we are interested in two lossy mutually coupled cavities, the Lindblad master equation, which couples a thermal bath with the dynamics induced by (1), takes the form: whereρ is the density matrix and γ is the loss rate. System (2) can be evolved numerically by computing quantum trajectories [6], where the wavefunction is evolved with a Monte Carlo algorithm. Such a computation is performed for a non-Hermitian effective Hamiltonian where the non-Hermitian terms account for nullmeasurement back-action; moreover, Eq. (3) is expressed in a co-moving frame with respect to the frequency ω p of the driving field, and Δ = ω c − ω p is the detuning.
The coherent evolution is interrupted by 'jumps', which correspond to photon emission from the cavities implemented by the action ofâ 1 orâ 2 on the wavefunction. The timing of these jumps is stochastic, and each trajectory corresponds to a particular unravelling of the master equation given by system (2). This approach is quite powerful as one can obtain temporal information regarding observables, which mimics the situation of an observer taking measurements of the outputs of the cavities; see [6] and Appendix A for more details. We consider here the observables a † 1 a 1 and a † 2 a 2 , the photon number expectation values in each of the two cavities.
In the co-moving frame with respect to the driving frequency ω p , the mean-field approximation is derived from system (2) by considering the expectations (mean-values) α 1 = â 1 and α 2 = â 2 under the assumption that their product factorises, meaning that, for example, â 1â2 = â 1 â 2 ; see [31] for details. It is mathematically convenient to rescale time and the variables of system (4) by introducing for the rescaled complex field electric amplitudes A and B, where System (5) is the semiclassical ODE model of the open Bose-Hubbard dimer, and we will investigate how its bifurcations manifest themselves in the evolution as described by the Lindblad master equation (2) for small photon numbers in the two cavities. Note that we are considering here the symmetric case, where both cavities are pumped with the same intensity F ; hence, system (2) is invariant under the exchange of the two cavities given by (â 1 ,â 2 ) → (â 2 ,â 1 ). This reflectional symmetry manifests itself as Z 2 -equivariance [19] of the semiclassical system (5) with respect to the map (A, B) → (B, A).
Mathematically, system (5) is a system of four realvalued ODEs, which can be written out either in terms of real and imaginary parts or amplitudes and phases of both A and B; hence, it defines a vector field with a four-dimensional phase space. Even though system (5) looks deceptively simple, the fact that its phase space is of dimension four implies that complex behaviour must be expected. Already continuous-time systems with three-dimensional phase spaces, such as the wellknown Lorenz [38] and Rössler systems [47], are known to be able to support chaotic dynamics and associated transitions and bifurcations; see, for example, the textbooks [21,33,50]. Specifically for system (5), a detailed study in [17] has uncovered a complex bifurcation structure with different types of chaotic dynamics in well characterized regions of phase space; these are organised by bifurcations involving homoclinic and heteroclinic connections of different kinds, some of which have not been studied before.
The sign ξ(U ) = sign(U ) appears in (5) to allow for convenient comparisons with different instances of the Bose-Hubbard dimer in the literature; however, the transformation implies that all results for ξ(U ) = 1 directly translate to those for ξ(U ) = −1. From now on, we set ξ(U ) = sign(U ) = 1 and consider positive κ as the experimentally relevant case [16,23,24] whose bifurcations have been studied in considerable detail in [17,18]. More specifically, we also set δ = −4.5 and consider the one-parameter bifurcation diagram in the pump strength f , for which one finds different types of phase transitions, including symmetry breaking, multistability and the onset of periodicity.
We compare the dynamics of the semiclassical ODE model (5) with quantum trajectory computations, requiring that we set (J, Δ, U, γ) = (−3.5, 4.5, 0.5, 2.0) in system (2) so that the two sets of parameter values agree. Regarding pump strength F , we note that any scaling of the form with μ > 0 does not change the value of f in system (5) as given by (6); hence, the rescaling (7) leaves the observed semiclassical behaviour unchanged. On the other hand, the introduction of the scaling parameter μ in the quantum optical system (2) appears in its effective Hamiltonian aŝ More precisely, by changing the scaling parameter μ, one obtains a quantum system operating with on-site energy U μ = U/μ and pump amplitude F μ = √ μF . In light of scaling (7), the μ-dependent effective Hamilto-nianĤ μ eff has the exact same semiclassical approxima- tion for any μ, namely, the semiclassical ODE model (5) derived from system (2). Nevertheless, in the quantum regime the scaling factor μ plays a vital role: varying it allows one to simulate a situation of a quantum system that supports different numbers of photons without changing the nature of its semiclassical limit. More specifically, the expected number of photons trapped in the cavities scales with μ with respect to the intensities predicted by the mean field model, namely, as a † 1 a 1 ≈ μ|A| 2 and a † 2 a 2 ≈ μ|B| 2 . In particular, increasing μ increases the number of photons and reduces the importance of cavity field fluctuations, meaning that the quantum system is closer to its semiclassical limit. In what follows, we investigate the predictive power of the semiclassical ODE model (5) by simulatingĤ μ eff in (8) for increasing values of μ; this is achieved by performing quantum trajectory simulations with the method detailed in Appendix A. We remark that this approach has been used in [37] to investigate system (1) for the case J > 0.

Bifurcation diagram and tracked quantum trajectory simulation
We first compare in Fig. 1 the one-parameter bifurcation diagram in the pump current f of system (5) for (δ, κ) = (−4.5, 3.5) with quantum trajectories of system (2) computed for slowly increasing pump current with corresponding parameter values. Figure 1(a) shows the branches and bifurcations of the equilibrium and periodic solutions of system (5) in terms of their intensities |A| 2 and |B| 2 over the relevant range of the pump strength f . Here, one finds stable equilibria along blue curves, while cyan and orange curves represent unstable equilibria with different numbers of positive eigenvalues. The green curves represent stable periodic solutions; specifically, these curves trace out the maxima (only) in |A| 2 and |B| 2 of these periodic solutions. Throughout, there exists a symmetric equilibrium with |A| 2 = |B| 2 and arg (A) = arg (B), and one finds a number of bifurcations that give rise to other branches of solutions. The f -range is divided in this way into intervals (i)-(vii) with different qualitative behaviour as follows: (i) f = 0 to the pitchfork bifurcation P 1 : there exists only the symmetric equilibrium and it is stable. (ii) P 1 to the Andronov-Hopf bifurcation H 1 : the symmetric equilibrium is now unstable and a pair of stable asymmetric equilibria exist; these bifurcate at P 1 with |A| 2 > |B| 2 and |A| 2 < |B| 2 , respectively. (iii) between the Andronov-Hopf bifurcations H 1 and H 2 : the asymmetric equilibria are unstable, and there is now a pair of stable periodic orbits that emerge and disappear at H 1 and H 2 -one near each equilibrium with |A| 2 > |B| 2 and |A| 2 < |B| 2 , respectively.
(iv) H 2 to the saddle-node bifurcation of asymmetric states S * 1 : the asymmetric equilibria are again stable and the situation is as in interval (ii).
(v) between the saddle-node bifurcations of asymmetric states S * 1 and S * 2 : two additional pairs of asymmetric equilibria exist, one stable and one unstable, so that there are now two pairs of stable asymmetric equilibria. The pair of unstable equilibria emerges or disappears with the respective pair of stable periodic orbits at the points S * 1 and S * 2 . (vi) S * 2 to the pitchfork bifurcation P 2 : there is again a single pair of stable asymmetric equilibria as well as the unstable symmetric equilibrium, as in intervals (ii) and (iv). (vii) beyond P 2 : the symmetric equilibrium is the only solution and stable again, as in interval (i).
shows the realisation of one quantum trajectory of system (2) as F μ with μ = 3.0 is increased linearly at a slow rate of approximately 0.216 per unit of time t, where (J, Δ, U, γ) = (−3.5, 4.5, 0.5, 2.0). The quantum trajectory is shown in terms of the two observables a † 1 a 1 and a † 2 a 2 , and superimposed is the bifurcation diagram from panel (a) after the corresponding scaling by μ of the intensities |A| 2 and |B| 2 ; here stable branches are dark grey and unstable branches light grey. Notice in Fig. 1(b) how both observables follow the branch of symmetric equilibria in interval (i) and then split into an asymmetric situation when the semiclassical system exhibits the pitchfork bifurcation P 1 . Near and beyond this transition, in intervals (i) and (iii), there are increased fluctuations that include switching between which of a † 1 a 1 or a † 2 a 2 is the larger observable. Past the Andronov-Hopf bifurcation H 2 , in interval (iv), the quantum trajectory is clearly localised at one of the asymmetric stable equilibria; more precisely, the one with a † 2 a 2 > a † 1 a 1 . As F μ increases, the quantum trajectory jumps to a different asymmetric equilibrium, namely, in between the asymmetric saddlenode bifurcations S * 1 and S * 2 , where two pairs of asymmetric equilibria exist, which is interval (v). As F μ increases further, a † 1 a 1 and a † 2 a 2 remain localised near this new asymmetric equilibrium in interval (vi). The two observables then come together near the pitchfork bifurcation P 2 and then remain localised at the symmetric equilibrium that is stable again in parameter interval (vii).
Overall, Fig. 1 illustrates that the slowly ramped quantum trajectory effectively follows the stable branches of the one-parameter bifurcation diagram, meaning that it remains localised near one of the stable equilibrium solutions, but with a considerable level of fluctuations. Increased sensitivity and jumps induced by these fluctuations are observed near bifurcation points, in the F μ -range between P 1 and H 2 as well as that bounded by S * 1 and S * 2 with multistability between different asymmetric equilibria.

Quantum trajectories at specific values of f
We proceed by investigating how these observed properties of the quantum trajectory manifest themselves in the parameter intervals with different limiting semiclassical dynamics; specifically, at the values of f = F μ / √ μ in intervals (i)-(vi), as indicated by the dashed vertical lines in Fig. 1(a). Here, we consider two values of the scaling factor μ, that is, two different photon numbers, namely, μ = 1.0 and μ = 3.0. More specifically, we present for each case: (1) temporal traces for a single realisation of a quantum trajectory of system (2) represented by the observ-ables a † 1 a 1 and a † 2 a 2 , shown over the corresponding range (that depends on μ) with the respective equilibria of system (5).
(2) the associated temporal traces of the ratio which is a measure of the validity of the factorisation property of the quantum system used to derive the semiclassical ODE; note that O = 1 means that the factorisation is exact.  (a1) and (b1) show the temporal trace of the observables â † 1â 1 (red curve) and â † 2â 2 (blue curve), and panels (a2) and (b2) show the associated evolution of the . Panels (a3) and (b3) show histograms for a 200 × 200 grid in the (|A| 2 , |B| 2 )-plane constructed from three different quantum trajectories; the data is symmetrised and the diagonal symmetry line is shown in grey. Also shown are the equilibria of system (5), as dashed lines in the temporal traces and as dots in the histograms, where colour reflects their stability as in Fig. 1 (a1) (b1) the respective equilibria and periodic solutions of system (5). The histograms are constructed from three different realisations of quantum trajectories by sampling 50,000 equidistant points in the time interval (0, 10 4 ); moreover, the invariance of system (2) under the permutation of the sites is used to double the number of points and to correct asymmetric bias introduced by the sampling.

Comparison in intervals (i) to (iii)
In interval (i), there exists only the single symmetric stable equilibrium of the semiclassical ODE (5), which is stable and attracts all initial conditions. As Fig. 2 shows, the respective quantum trajectory stays close to this equilibrium but is subject to clear fluctuations. This localisation is illustrated in panels (a1) and (b1) with the temporal traces of the two observables â † 1â 1 and â † 2â 2 , which can be seen to fluctuate around the corresponding equilibrium intensity-value. Here the ranges have been chosen to agree with the scal-ing by μ so that a direct comparison of the observables can be made, including with their limiting semiclassical behaviour. Notice that the fluctuations are relatively larger for μ = 1 in panel (a1) than those for μ = 3 in panel (b1). Similarly, the observable O, while subject to fluctuations in both cases, is on average further from its limiting value of 1.0 for μ = 1 in panel (a2) compared to μ = 3 in panel (b2). It is interesting to note that, despite the factorisation assumption being relatively inaccurate in this regime, the two-dimensional histograms are still well centred around the stable equilibrium. The twodimensional histograms in the respective ranges of the ( â † 1â 1 , â † 2â 2 )-plane in Fig. 2(a3) and (b3) show distributions that are well centered around the stable equilibrium on the symmetry axis. Note that the colour map is scaled to account for the dependence of the size of the bins on μ; this also allow for a direct comparison of the (relative) heights of the histograms for different values of μ. There is considerable spread due to the fluctuations in the quantum trajectories, which are indeed comparable with the temporal traces in panels (a1) and (b1). The spread is smaller and more symmetrical around the equilibrium for μ = 3 in panel (b3) compared to that for μ = 1 in panel (a3). This illustrates that the distribution becomes more Gaussian with a smaller variance as μ is increased. In interval (ii), past the pitchfork bifurcation P 1 , the symmetric equilibrium of system (5) is now unstable, and there is a pair of stable asymmetric equilibria, each with their own basin of attraction. Figure 3 shows that quantum trajectories of system (2) display switching between the two asymmetric states in between (quite short) epochs of localisation near one of them. The switching appears to be dominant in the temporal trace for μ = 1 in panel (a1), and epochs of localisation (while still short) are visible more clearly for μ = 3 in panel (b1). Overall, the role of fluctuations seems to be much more important here than in interval (i), as they drive switching between the two stable solutions. These observations are represented in the twodimensional histograms in panels (a3) and (b3) by the fact that the distributions are now bimodal, quite broad and not sharply focused around the two stable asymmetric equilibria. Notice also the existence of a 'bridge' between the areas of localisation near the asymmetric equilibria, which reflects the likely route for the switching driven by the quantum fluctuations. Again, the histogram is less broad, and its features are crisper for μ = 3 in panel (b3) compared to μ = 1 in panel (a3). Figure 4 illustrates the situation in interval (iii), in between the two Hopf bifurcation points H 1 and H 2 , with the new feature of a pair of attracting periodic solutions of system (5) near the now unstable asymmetric equilibria. The quantum trajectories of system (2) still display switching between these two asymmetric periodic states with epochs near one of them, with a considerable level of fluctuations. As before, the level of fluctuations is relatively higher for μ = 1 in panels (a1) and (a2) than for μ = 3 in panels (b1) and (b2). Comparison with Fig. 3 shows that the fluctuations are larger compared to the situation in interval (ii). In particular, the fluctuations during epochs of localisation are now larger, since the quantum trajectories are √ μ = f = 6.0 as represented in Fig. 2; also shown in the histograms is the pair of stable asymmetric periodic solutions (green closed curves) no longer attracted to a steady state. The histograms in Fig. 4(a3) and (b3) show distributions that illustrate the properties of quantum trajectories differently. While the maxima are near the stable periodic solutions of the semiclassical ODE, there appears to be no definite fingerprint of periodicity of the quantum trajectories in the ( â † 1â 1 , â † 2â 2 )-plane. Notice that there is again a clear 'bridge' of preferred switching, which is considerably sharper for μ = 1 compared to μ = 3, but the histograms do not appear to identify the pair of periodic attractors.
The above discussion shows that the frequent switching between the two localised oscillations obscures the possible periodicity of the observables â † 1â 1 and â † 2â 2 . To identify the fingerprint of these semiclassical periodic oscillations in the quantum realm, we now show in Fig. 5 temporal traces of the observable S o = â † 1â 1 + â † 2â 2 of the quantum trajectories as well as their spectra, for μ = 1, μ = 3 and also for μ = 20. The observable S o (subject to the same scaling by μ) is the quantum analogue of the total intensity S = |A| 2 + |B| 2 of the semiclassical ODE. Due to the symmetry properties of the Bose-Hubbard dimer, switching between symmetric states manifests itself much less in S o . More specifically, due to its invariance under the permutation of the two sites, S o minimises fluctuation transients driven by the quantum system. This is why one can observe signs of periodicity in the temporal traces in Fig. 5(a1)-(c1). The temporal trace of the semiclassical periodic orbit is shown in Fig. 5(d) for comparison; note that it is close to being sinusoidal, which is due to the periodic orbit still being close to the Hopf bifurcation. As μ is increased, fluctuations are reduced, and the periodicity in the temporal trace of S o becomes crisper. Indeed, the temporal trace for μ = 20 in panel (c1) is recognised as a 'noisy version' of the periodic signal in panel (d) and, hence, clearly contains fingerprints of the semiclassical periodic solution. This observation is quantified by power spectra |F(S o )| of the respective temporal traces. Already for μ = 1 in panel (a2) the spectrum shows a recognis- √ μ = f = 6.0.

Panels (a1), (b1) and (c1)
show temporal traces of the observable So = â † 1â 1 + â † 2â 2 of quantum trajectories of system (2) with μ = 1, μ = 3 and μ = 20, respectively. Panels (a2), (b2) and (c2) show the corresponding power spectra |F (So)| (green data); here the vertical purple line indicates the frequency of the corresponding periodic solution of system (2) able peak near the main frequency of the semiclassical periodic temporal trace. For μ = 3 in panel (b2), the spectrum is sharper and its frequency closer to that of the semiclassical oscillation, and this is even more the case for μ = 20 in panel (c2). We conclude that Fig. 5 clearly shows the emergence of periodicity in the quantum realm provided μ, that is, the photon number, is taken to be sufficiently large (but still moderate). We remark that this phenomenon has also been observed recently in quantum trajectory simulations of the unbalanced Dicke model [49].

Comparison in intervals (iv) to (vi)
In interval (iv), past the second Hopf bifurcation H 2 , the two asymmetric equilibria are again stable and the only attractors of system (5). In other words, the situation is topologically the same as in interval (ii). However, as Fig. 6 shows, we find marked differences in the observed behaviour of quantum trajectories of system (2). The temporal trace for μ = 1 in panel (a1) show long epochs, where the quantum trajectory is localised near one of the stable steady states, with much more occasional switchings between them. For μ = 3 in panel (b1), the relative strength of quantum fluctuations is now so low that not a single switching occurs in the time window presented. Comparison with Fig. 3 shows that the overall level of fluctuation is comparable in intervals (ii) and (iv). However, the observable O is now seen to fluctuate around unity in panels (a2) and (b2), meaning that the factorization assumption is reasonable in interval (iv). Notice further that switching events manifest themselves in the observable O in panel (a2) as sudden larger spikes away from its average. The clear observation that quantum trajectories linger much longer near one of the two stable equilibria in interval (iv) is explained by the fact that these equilibria are more attracting and also further apart from each other and from the unstable symmetric equilibrium from which they bifurcated. More specifically, the two stable equilibria represent a situation of extreme symmetry breaking, where one site has practically all photons of the overall coupled system, while the other has near-zero photons-this in spite of the fact that both sites are pumped symmetrically. This distance between the two attractors and the increased residence times of quantum trajectories are illustrated very clearly by the histograms in panels (a3) and (b3). Already for μ = 1 in panel (a3), the distribution is very bimodal with pronounced peaks near the two attracting equilibria system (5), which are very close to the coordinate axes owing to the fact that one of the two intensities is practically zero; the 'bridge' corresponding to preferred switching is now much weakened. For μ = 3 in panel (b3), there is no longer a discernible bridge due to the occurrence of very few transitions between the two attractors. The histogram is now quite sharply focused on the pair of stable asymmetric equilibria, with very well defined peaks. That notwithstanding, there are still sufficiently many switchings due to quantum fluctuations to ensure that the computed histogram captures both attractors. The behaviour in interval (v), in between the two saddle-node bifurcations S * 1 and S * 2 , is characterised by the existence of an additional pair of stable asymmet-ric equilibria of system (5), as well as a pair of unstable asymmetric equilibria. This means that there are now a total of four attractors with their respective basins of attraction. The additional stable equilibria are characterised by a relatively small imbalance between the sites compared to the other pair, which still represent the scenario where practically all photons are at one of the two sites. Figure 7 shows that quantum trajectories of system (2) remain localised near one of the four stable equilibria for a certain amount of time and then switch to being localised near another stable equilibrium and so on. The temporal trace for μ = 1 in panel (a1) shows quite a number of switchings. While for μ = 3 in panel (b1), residence times are longer, and the number of switchings per time interval is decreased. Interestingly, the overall level of fluctuation is relatively quite low compared to earlier cases. Note also that, as seen in panels (a2) and (b2), localisation near the new equilibria with smaller values of the intensity is associated with especially low fluctuations, namely around an average value of O = 1.0, implying the factorization assumption is accurate in this regime. These attractors feature significant photon numbers in both cavities, such that the role of quantum fluctuations is less important. However, while the other attractors feature a large population in one cavity, the other, almost empty, cavity is strongly impacted by fluctuations. The four different attractors are also clear features of the histograms in panels (a3) and (b3). The distributions show clear peaks near each of the stable equilibria, which are considerably crisper for μ = 3 compared to μ = 1. Moreover, the histograms are considerably larger near the pair of equilibria, where one of the intensities is practically zero. Note that switching events follow 'weak bridges' between the attractors, meaning that switching between neighbouring attractors in the ( â † 1â 1 , â † 2â 2 )plane are by far the ones that are most likely to occur. Figure 8 shows the situation in interval (vi), where there are again only two stable asymmetric equilibria of system (5), namely, the ones near the symmetry line with a relatively small imbalance between the sites. As panels (a1) and (b1) show, the quantum trajecto-ries of system (2) remain near both of these equilibria, with negligible residence times near each attractor and many switchings per time interval. Residence times are notably larger for μ = 3 than for μ = 1, but remain very small for either case compared to the (topologically equivalent) situation in interval (iv) in Fig. 6. The level of fluctuations appears to be relatively small, and the observable O in Fig. 8(a2) and (b2) remains very close to 1.0 throughout; again, any noticeable fluctuations of O appear to be due to switchings between the two stable equilibria. However, the role of these fluctuations is important, as the frequent switching means that the associated histograms in panels (a3) and (b3) are characterised by a single peak with a maximum at a point in between the two stable asymmetric equilibria. In other words, the two nearby attractors are not resolved: while the histogram for μ = 3 in panel (b3) appears to be more elongated around the two attracting equilibria, it still has only a single peak. This is a direct reflection of the very low residence times of localisation as measured by the observables â † 1â 1 and â † 2â 2 . As we will see below, it is possible to detect the (weak) localisation of the quantum trajectories already for μ = 3, namely, by considering the difference between these two observables; see already Fig. 9(b2). Figure 9 provides an overview of how the statistical properties of quantum trajectories of system (2) develop with the pump strength F μ , and how this compares with the corresponding bifurcation diagram of the semiclassical ODE (5). Here, we again consider the two cases μ = 1 and μ = 3. In this representation, histograms are shown for F μ starting from 2.0μ in steps of 0.5μ and up to 15.75μ, where the maximum of each histogram is scaled to the F μ -size of 0.5μ. As for the two-dimensional histograms shown in previous figures, these histograms are constructed from three different realisations of quantum trajectories by sampling 50,000 equidistant points in the time interval (0, 10 4 ), where we now use 80 uniform bins over the shown range of the respective observable. The invariance of system (2) under the permutation of the sites is again used to double the number of points and symmetrise the data; the maxima are then scaled to 0.5μ. Figure 9 shows histogram plots for two different observables. Panels (a2) and (b2) show histograms for â † iâ i , which corresponds to the projection of the respective two-dimensional histogram onto the â † 1â 1axis (or equivalently the â † 2â 2 -axis). To obtain additional information regarding semiclassical fingerprints in the distributions, Fig. 9(a2) and (b2) show histograms for the difference D o = â † 1â 1 − â † 2â 2 ; note that considering this observable corresponds to the projection of the respective two-dimensional histogram onto the antidiagonal and implies the symmetry of panels (a2) and (b2) with respect to the F μ -axis.

Evolution of histograms with F µ
The evolution of the histograms for μ = 1 is shown in Fig. 9(a1) and (a2). There are noticeable changes in the statistical properties of the two observables associated with the transitions through the different bifurcations, while the distributions remain more or less the same in the intervals (i) to (vi) that are covered by the shown range of F μ ; compare with Fig. 1. Starting at low F μ in interval (i), one first observes that the histogram widens as the first pitchfork bifurcation is approached and subsequently becomes bimodal in interval (ii), with a pair of peaks near the two stable asymmetric equilibria. Notice the large component of the distribution in between these two equilibria, which corresponds to frequent switchings between them. This 'bridge' and, hence, the number of switchings in the time interval clearly become much less pronounced past interval (iii). The existence of stable periodic orbits, on the other hand, is not evident in the histograms. When the second Hopf bifurcation is reached, the distribution is strongly bimodal and remains clearly localised near the re-stabilised asymmetric equilibria throughout interval (iv). Notice that in Fig. 9(a1) for the observable â † iâ i , the solution is less well resolved near the upper branch compared to the lower branch due to the larger fluctuations for the site with more photons; this issue clearly does not arise for the symmetric observable D o in panel (a2). The strong localisation extends well into interval (v) with the additional pair of stable equilibria. There is a somewhat gradual change of the histogram to localisation around the other pair of stable equilibria in this interval. Notice that the distinction between these two equilibria is quite weak initially and quickly becomes nonexistent, with a distribution with a maximum in between the two attractors, even for the observable D o .
The evolution of the histograms for μ = 3 in Fig. 9(b1) and (b2) is quite similar, but the distributions are more clearly resolved, that is, more concentrated at the respective attractor. However, there are noticeable differences in the upper range of F μ . In interval (v), the switching to a different pair of attractors now manifests itself as a quite sudden change of the histogram; this reflects the scarcity of switchings in a finite time series. Moreover, the two new asymmetric equilibria are now distinguished by the histogram, especially clearly for the observable D o in panel (b2). Overall, we conclude that the histogram plots of Fig. 9, especially those in panels (a2) and (b2) for the symmetric observable D o , provide a good summary of how the statistical properties of quantum trajectories change both with the pump strength F μ as well as with increasing numbers of photons as represented by the scaling parameter μ. Indeed, this representation agrees with the results in Sect. 4 regarding the behaviour for representative values of F μ in the intervals (i) to (vi)-but it also provides insight into how the distributions of quantum trajectories change from interval to interval, as semiclassical bifurcations are encountered.

Antibunching and entanglement
We now investigate whether there is evidence of quantum phenomena in system (2) for low photon numbers as the pump strength F μ is varied. To this end, Fig. 10 shows the two indicator functions min{g (2) 11 (0), g (2) 22 (0)} and E over the relevant range of F μ , where the semiclassical bifurcations are shown and the associated intervals (i) to (vii) of different behaviour are highlighted.
Antibunching of photons refers to the emission of a photon reducing the probability of a subsequent emission, and is a strictly quantum phenomena. In the present context, it can be detected by the condition that min{g (2) 11 (0), g (2) 22 (0)} < 1, where the function g (2) ii (0) is the second-order correlation function for the light in cavity i g (2) ii Plotting the smaller of the time averaged secondorder correlation functions of the two sites, g 11 (0) and g (2) 22 (0), for a single quantum trajectory allows us to identify antibunching [39]. Figure 10(a) shows timeaveraged values of min{g (2) 11 (0), g (2) 22 (0)} for μ = 0.5, μ = 1.0, μ = 2.0, and μ = 3.0 at different values of F μ for single trajectories of length t = 500. In interval (i), both sites display slight bunching, less so for larger μ, implying a thermal photon number distribution. Past the pitchfork bifurcation P 1 and the appearance of asymmetric attracting states of system (5), we find min{g (2) 11 (0), g (2) 22 (0)} < 1: consistently in intervals (iii) to (vii) for all shown values of μ, and including interval (ii) for μ = 0.5 and μ = 1.0. Hence, the site with more photons is always anti-bunched from about F μ = 5.0 and higher. As one would expect, antibunching is reduced as μ increases and the closer the system is to its semiclassical limit, and min{g (2) 11 (0), g (2) 22 (0)} appears to approach the value 1.0 as μ is increased.
To investigate the existence of entanglement between the sites, we consider the time-averaged von Neumann entropy of one of the cavities, given by [40,43] where ρ 1 (t) is the density matrix for site 1 at time t (or equivalently site 2) and the overline represents time averaging. Figure 10(b) shows E as a function of F μ for μ = 0.5, μ = 1.0, μ = 2.0, and μ = 3.0, as computed from the same quantum trajectories described above. Entanglement is identified by the condition that E > 0, which means that the states are not 'pure', and we conclude that there is entanglement throughout the entire range of F μ . Notice that the level of entanglement as measured by E increases in interval (i) and then, past the first pitchfork bifurcation P 1 , reaches and stays on a plateau throughout interval (ii) to (v). The indicator E then increases quite steeply in interval (vi) with maximal multistability between asymmetric attractors, reaches a maximum near the second pitchfork bifurcation P 2 and then decreases equally steeply in interval (vii) where the symmetric equilibrium is again the only attractor of the limiting system (5). The von Neumann entropy in the plateau appears to show slowly reducing entanglement as μ increases, suggesting that entanglement goes away in the thermodynamic limit. However, the maximum of E appears to become larger and more narrow with increasing μ. It is a known phenomenon that, near certain types of phase transitions, there is a divergence in entanglement at the critical point in the thermodynamic limit [35,36,41,45]. Whether this is the explanation for the sharpening of the von Neumann energy near P 2 remains an interesting question for future research.

Conclusions and outlook
The case study of the open two-site Bose-Hubbard dimer presented here shows that it is possible to identify recognisable fingerprints of phase transitions-that is, bifurcations-of the limiting semiclassical (mean-field) model in the quantum realm. More specifically, we considered the case of negative intermode coupling when the semiclassical model features a transition, as the pump strength is increased, from symmetric dynamics via symmetry breaking at a pitchfork bifurcation to oscillatory dynamics and then to multistability between different types of asymmetric states. These features were recognised reliably in the statistical properties of different observables of quantum trajectories, even quite far from the semiclassical limit, that is, for low numbers of photons at each site.
These theoretical results are fundamental in nature and suggest that it may be possible to find quantum signatures of such rather complex nonlinear phenomena of the Bose-Hubbard dimer in different experimental systems, including condensates in semiconductor microcavities [1,34,46], superconducting circuits [13,44] and photonic crystal nanocavities [16,18,23,24]. In semiconductor microcavity condensates, oscillatory behaviour [1,34] and multi-stability under asymmetric driving [46] have been observed. Moreover, superconducting circuits have been built that act as tuneable degenerate and non-degenerate amplifiers under asymmetric drive [13]. Finally, for two coupled photonic crystal nanocavities [23,24], symmetric pumping has recently been achieved reliably up to a moderate level of the pump strength, enabling the experimental verification of spontaneous symmetry breaking [16] in good agreement with a bifurcation study of the semiclassical model [18]. We, therefore, believe that more complicated dynamics with additional levels of multistability-including different types of localised and non-localised chaotic dynamics [17]-may well be identifiable experimentally in these systems; especially, in photonic crystal structures due to their small size and low-photon operation. Such experiments for the Bose-Hubbard dimer, or other quantum systems such as the Dicke model [2,3,25,30,48,49,53], are challenging but provide opportunities for studying how semiclassical chaos arises in quantum systems [4,15].
Finally, we remark that system (5) is part of a bigger class of Z 2 -equivariant equations that describes two identical optical systems that are being run and coupled symmetrically. This class contains ODE models, such as system (5) studied here, but also models in the form of partial differential equations (PDEs). As an example of the latter, we mention the system of two coupled Lugiato-Lefever equations (LLEs) that describe two counter-propagating modes in a micro-ring resonator. Recent results in [52] have shown good agreement between experiments and the predictions of the coupled LLEs for the case when the spatial component can be neglected; the PDE model then reduces to an ordinary differential equation for the complex electric fields of the two counter-propagating modes, which, at a first glance, looks very similar to system (5). However, the coupling between these two modes arises from the Kerr nonlinearity of the micro-ring resonator and is, therefore, nonlinear [52]-unlike the coupling in the Bose-Hubbard dimer which is linear [23,24]. This difference in the nature of the coupling is responsible for substantial differences in the organisation of observable behaviour. The study of emerging phenomena that are exclusive to the case of nonlinear coupling, including their quantum optical description, is a promising direction of ongoing work. A related interesting challenge for future research is the study of the general case of Bose-Hubbard systems with multiple sites [12] and different types of coupling.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.

Appendix: The quantum trajectory method
Quantum trajectories offer a method for calculating the evolution of a dissipative quantum system [6]. The main idea is that stochastically generated trajectories mimic single experimental runs, where energy dissipated into the environment is measured and recorded rather than lost. Each trajectory thus constitutes a particular unravelling of the master equation, such that an ensemble average over many trajectories reproduces the master equation result.
Quantum trajectories recast each time step in a master equation as one of two options: discrete quantum jumps and coherent non-unitary evolution. The former corresponds to the measurement of a dissipated photon in the environment and the latter to a null measurement. Non-unitary evolution is performed with an effective Hamiltonian where the additional term arises from the backaction of measuring zero dissipated photons. Discrete jumps are performed by the action of the jump operators {Ôj}. Evolution from an initial pure state |ψ(0) is calculated in time steps δt, where at each time step, |ψ(t+δt) is deter- MCstep 3A: If rn < pn then the time step is designated as a quantum jump. Which operator is used to perform the jump is decided stochastically, with the action ofÔj happening with probability: .
The state is updated as .
MCstep 3B: If rn > pn then the time step is designated as coherent evolution. The state is updated as .
We employ this method to produce long quantum simulations of the effective μ-dependent HamiltonianĤ μ eff of the open Bose-Hubbard dimer, as given by Eq. (8). The resulting statistical properties can then be compared with the prediction from the bifurcation analysis of the semiclassical ODE model (5). For these quantum simulations, the Hamiltonian, with scaling parameter μ, is, therefore, given byĤ jâ jâj + √ μF âj +â † j and the jump operators areâ1 andâ2, with rates γ1 = γ2 = γ, such that coherent evolution is performed with the non-Hermitian effective Hamiltonian whereÔ1 =â † 1 andÔ2 =â † 2 are our jump operators; compare with (9).