Driven Bose-Hubbard dimer under nonlocal dissipation: A bistable time crystal

TL;DR

This paper explores how nonlocal dissipation induces multistability and time-crystalline behavior in a driven Bose-Hubbard dimer, revealing novel critical phenomena and symmetry-breaking in open quantum systems.

Contribution

It demonstrates the emergence of a period-bistable time crystal and multistability due to nonlocal dissipation, combining semiclassical and quantum analyses.

Findings

Identification of amplitude bistability and symmetry breaking.

Formation of a period-bistable time crystal with Josephson-like oscillations.

Spectral analysis of the Liouvillian reveals critical behavior in the open quantum system.

Abstract

We investigate the critical behavior of the open coherently-driven Bose-Hubbard dimer under nonlocal dissipation. A conserved quantity arises from the nonlocal nature of the dissipation, rendering the dimer multistable. In the weak-coupling semiclassical limit, the displayed criticality takes the form of amplitude bistability and breaking of spatial and temporal symmetries. A period-bistable time crystal is formed, consisting of Josephson-like oscillations. Mean-field dynamics and quantum trajectories complement the spectral analysis of the Liouvillian in the approach to the semiclassical limit.

Equations

$${\partial }_t\hat{\rho }=-i[\hat{\mathcal{H}},\hat{\rho }]+\gamma \mathcal{D}[{\hat{a}}_1+{\hat{a}}_2](\hat{\rho })\equiv \mathcal{L}(\hat{\rho }).$$

$$\begin{array}{rl}\hat{\mathcal{H}}=\sum _{i=1,2}-\Delta & {\hat{a}}_i^{†}{\hat{a}}_i+U{\hat{a}}_i^{†}{\hat{a}}_i^{†}{\hat{a}}_i{\hat{a}}_i+F({\hat{a}}_i^{†}+{\hat{a}}_i)\\ & -J({\hat{a}}_1^{†}{\hat{a}}_2+{\hat{a}}_1{\hat{a}}_2^{†}),\end{array}$$

$${\hat{\rho }}_{\text{ss}}(t)={\hat{r}}_{0,1}+\sum _{d=2}^{D_0}{c}_{0,d}{\hat{r}}_{0,d}+\sum _{d=1}^{D_1}(c_{1,d}{e}^{i|{\lambda }_1|t}{\hat{r}}_{1,d}+\text{H.c.}),$$

Full text

Bose-Hubbard dimer, dissipative quantum phase transitions, nonlocal dissipation, multistability

I Introduction

Phase transitions substantiated by a mean-field bifurcation have been intimately tied to the evolution of quantum optics. Characteristic cases include the Dicke phase transition Kirton et al. (2019), co-operative resonance fluorescence Narducci et al. (1978), and the laser DeGiorgio and Scully (1970). The study of quantum phase transitions in zero dimensions Hwang et al. (2015) far from thermal equilibrium reappraises the definition of the thermodynamic limit alongside the distinct role of quantum fluctuations, taking the form of an ostensible deviation from the mean-field dynamics Savage and Carmichael (1988); Carmichael (2015). Furthermore, the work of Ref. Kessler et al. (2012) has already drawn the analogy between quantum phase transitions in closed systems and dissipative phase transitions (DPTs) in open quantum systems, through analyzing the spectrum of the Liouvillian, while the properties of the steady state have been assessed with respect to the particular critical eigenvalue defining the Liouvillian spectral gap Minganti et al. (2018). To date, DPTs attract a significant and ongoing interest among the communities of cold atoms, circuit QED, and semiconductors (see, e.g., Refs. Dimer et al. (2007); Baumann et al. (2010, 2011); Brennecke et al. (2013); Fitzpatrick et al. (2017); Rodriguez et al. (2017); Fink et al. (2017, 2018)).

In that background, the interplay between coherence and dissipation renders the Bose-Hubbard dimer (BHD) a suitable candidate for assessing the role of conserved quantities and the associated symmetries in dissipative quantum phase transitions. A transition from intermittent entanglement to a persistently entangled state was demonstrated a few years ago for a dissipative BHD Pudlik et al. (2013). The standard phenomenological approach to dissipation amounts typically to the definition of a Liouville superoperator governing the system response, assessed against its coherent evolution. Statistical independence of the reservoirs the system is coupled to gives rise to the configuration of local dissipation, while coupling to one common bosonic reservoir generates nonlocal dissipation.

Critical behavior in the presence of a coherent drive, leading to symmetry breaking, has been very recently taken up for local dissipation when characterizing the quantum correlation properties of the BHD Cao et al. (2016); Casteels and Ciuti (2017). In many cases, local dissipation is well justified and often greatly simplifies the analysis of the transient and the steady state. In reality, however, different collective modes have different dissipation rates, giving rise to some amount of nonlocal dissipation.

In this work, we demonstrate the importance of the nonlocality of the dissipation in the open coherently-driven BHD. A continuous swapping symmetry together with a conserved quantity arises, leading to the occurrence of quantum multistability Albert and Jiang (2014). Alongside the familiar first-order DPT Casteels et al. (2017) associated with semiclassical amplitude bistability, with origins in the quantum Duffing oscillator Drummond and Walls (1980), the BHD breaks the spatial swapping symmetry as well as the time invariance of the steady state. By studying the closure of the Liouvillian gap as we approach a weak-coupling semiclassical limit, we show that a dissipative time-crystal Sacha and Zakrzewski (2018); Iemini et al. (2018); Wang et al. (2018); O’Sullivan et al. (2018); Tucker et al. (2018); Gong et al. (2018); Gambetta et al. (2019) with Josephson-like oscillations is formed. We find that two different time-crystalline periods coexist in the parameter regime of semiclassical bistability. We employ and compare the semiclassical Gross-Pitaevskii equations, quantum trajectories, and numerical solutions of the Master Equation in a truncated Hilbert space in order to fully appreciate the novel phases.

II The model

We consider a BHD whose dynamical evolution is governed by the Lindblad equation (with $\hbar=1$)

[TABLE]

Here, $\mathcal{D}[\hat{L}](\bullet)=\hat{L}\bullet\hat{L}^{\dagger}-(1/2)\{\hat{L}^{\dagger}\hat{L},\bullet\}$ is the dissipator, while the system Hamiltonian under coherent drive, in a frame rotating with the pump frequency $\omega_{p}$, reads

[TABLE]

where $\hat{a}_{i}$ is the $i$-mode bosonic annihilation operator, $\Delta=\omega_{p}-\omega_{c}$ is the detuning between the pump and the resonant mode frequencies, $\omega_{p}$ and $\omega_{c}$ respectively, $U$ is the interaction strength, $F$ is the driving amplitude, and $J$ is the intermode coupling. For nonlocal dissipation, $\hat{L}=\hat{a}_{1}+\hat{a}_{2}$. The driving term originates from the interaction of the two bosonic modes $\hat{a}_{1,2}$ with a coherent source treated as a (semi)classical field. Drive and decay only affect the bonding mode $(\hat{a}_{1}+\hat{a}_{2})/\sqrt{2}$, which couples to the antibonding mode $(\hat{a}_{1}-\hat{a}_{2})/\sqrt{2}$ solely via the nonlinearity. We rescale the pump amplitude and the interaction with the help of a generic parameter $N$ (the analog of the laser saturation photon number, for instance), as $F=\sqrt{N}\tilde{F}$ and $U=\tilde{U}/N$, such that $\sqrt{U}F$ is constant for any $N$. This allows us to consider a well defined weak-coupling thermodynamic limit, $N\to\infty$, where the boson number diverges ($\langle\hat{a}_{i}^{\dagger}\hat{a}_{i}\rangle\propto N$). In this limit, quantum fluctuations are negligible, and a semiclassical analysis is adequate Casteels and Ciuti (2017). Our system exhibits three important properties delineated below.

The first one is the presence of a continuous symmetry $\mathcal{L}(\hat{U}(\phi)\bullet\hat{U}(\phi)^{\dagger})=\hat{U}(\phi)\mathcal{L}(\bullet)\hat{U}(\phi)^{\dagger}$ generated by the unitary operator $\hat{U}(\phi)=e^{i\phi\hat{Z}_{2}}$ (for real $\phi$), with the swapping operator $\hat{Z}_{2}=\sum_{n_{1},n_{2}}\ket{n_{1},n_{2}}\bra{n_{2},n_{1}}$ (written in the Fock-state basis of the two individual modes). The swapping operator is a conserved quantity. This is only possible due to the nonlocality of the dissipation Nigro (2019) and is the origin of quantum multistability Albert and Jiang (2014). Since there are infinitely-many possible initial expectation values of $\hat{Z}_{2}$, there are also infinitely-many steady states. This is true for any value of $N$. For a finite $N$, the steady state reads $\hat{\rho}_{\text{ss}}=\hat{r}_{0,1}+c_{0,2}\hat{r}_{0,2}$, where $\hat{r}_{0,1}$ is a density matrix while $\hat{r}_{0,2}$ is not, and the coefficient $c_{0,2}=\operatorname{Tr}[\hat{Z}_{2}\hat{\rho}(0)]$ depends on the initial condition. This steady state is symmetric, i.e., $\hat{U}(\phi)\hat{\rho}_{\text{ss}}\hat{U}^{\dagger}(\phi)=\hat{\rho}_{\text{ss}}$, irrespective of the initial condition.

Secondly, the system we consider presents limit cycles 111Note that the limit cycles are not a consequence of having a nontrivial conserved operator ($\hat{Z}_{2}$). The same BHD model, but with antisymmetric drive, also exhibits limit cycles and only $\hat{\mathbb{1}}$ is conserved. in the limit $N\to\infty$. These are persistent periodic oscillations in the infinite time limit of the dynamical evolution. Recently, an open quantum system with this property has been termed a boundary time crystal Iemini et al. (2018). It is there argued that while the total Hamiltonian (of the boundary, the bulk, and their coupling interaction) is time-translation invariant, the boundary system presents limit cycles, thus breaking the global time-translation symmetry and forming time crystal with a nonrigid period that is allowed to change continuously as a function of the system parameters. Our system satisfies the same conditions. The dependence on the period of the coherent drive can be eliminated by a gauge transformation. This means that any emergent periodic response stems from a continuous breaking of time-translation invariance in contrast to the discrete fashion due to the presence of a Floquet-map eigenvalue in the unit circle Wang et al. (2018); Gong et al. (2018). The steady state of the BHD, which spontaneously breaks the time-translation invariance of the Liouvillian, assumes the form (for $N\to\infty$)

[TABLE]

where $\mathcal{L}(\hat{r}_{n,d})=\lambda_{n}\hat{r}_{n,d}$, $c_{n,d}$ are coefficients depending on the initial system density matrix $\hat{\rho}(0)$, and $\lambda_{1}$ is a purely imaginary eigenvalue responsible for the formation of a limit cycle 222See the Supplemental Material for further details, which includes Refs. Buca et al. (2019); Aleiner et al. (2012); Rayanov et al. (2015); Macieszczak et al. (2016); Plenio and Knight (1998); Breuer and Petruccione (2002); Hofer et al. (2017); Ciccarello (2018); Choi et al. (2017).. The upper limits $D_{0}$ and $D_{1}$ in the sums of the right-hand side in Eq. (3) are unknown. The index $d=1,...,D_{n}$ labels the degeneracy of the eigenvalue $\lambda_{n}$. In principle, one could prepare an initial condition such that the coefficients $c_{1,d}$ vanish and the system reaches a time-independent steady state. For a system prepared in such a manner, an infinitesimal disturbance would lead to $c_{1,d}\neq 0$ with some degeneracy $d$, causing the response to oscillate ceaselessly and forming a time crystal. This could in principle be observed in any time-delayed two-point measurement.

Finally, the system response is determined by semiclassical complex-amplitude bistability as $N\to\infty$. This is known to occur in the BHD under local dissipation for the same driving configuration as the one we consider here Casteels and Wouters (2017). Semiclassical bistability manifests itself as a region in the parameter space where two semiclassical fixed points, with different complex amplitudes, appear. In our model, we will show that two attractors, specifically limit cycles with different periods, coexist in this region.

III Results

In order to assess the system behavior as $N\to\infty$ we resort to the Gross-Pitaevskii (semiclassical) approach, in which the operators in the Heisenberg picture are replaced by their expectation values $\alpha_{i}=\langle\hat{a}_{i}\rangle$. We remark that our continuous swapping symmetry is intrinsically quantum and cannot be captured by the semiclassical picture, where it shows up merely as a discrete swapping symmetry $\alpha_{1}\leftrightarrow\alpha_{2}$. We compare the semiclassical predictions with the calculations in the quantum regime for increasing values of $N$.

In Fig. 1(a) we depict the steady state expectation value of the rescaled total number of bosons $\hat{n}_{\text{tot}}/N=(\hat{a}_{1}^{\dagger}\hat{a}_{1}+\hat{a}_{2}^{\dagger}\hat{a}_{2})/N$ as a function of the rescaled pump amplitude $\tilde{F}/\gamma=F/(\sqrt{N}\gamma)$. We present in color the quantum results for different values of $N$. The continuous (dashed) black lines indicate stable (unstable) semiclassical fixed points. One can see that as $N$ increases, the occupation number in the quantum regime approaches one of the two semiclassical stable branches. The transition between the two branches becomes increasingly sharper, suggesting the formation of a discontinuous jump, as one would expect from a first-order DPT Casteels et al. (2017); Minganti et al. (2018).

The basin of attraction of the stable fixed points is just the two-dimensional plane given by $\alpha_{1}=\alpha_{2}$ in the four-dimensional phase space defined by the real and imaginary parts of $\alpha_{1}$ and $\alpha_{2}$. A randomly chosen initial condition generates a trajectory that most certainly avoids the plane, converging to one of the families of limit cycles which revolve around the stable fixed points. The amplitude of the attained limit cycle is determined by the initial difference in the coherent states $\alpha_{1}(0)-\alpha_{2}(0)$, although its period is solely dependent on the set of parameters entering in Eq. (1), except in the region of bistability, where limit cycles with two different frequencies have different basins of attraction. Interestingly, the limit cycles break the swapping symmetry. They are the so-called Josephson oscillations. A typical example of these oscillations is depicted in Fig. 1(b).

In Fig. 1(d), we show the five eigenvalues whose real parts are closer to zero (defining the sequence of the gaps), for increasing values of $N$. We distinguish them using different colors and markers according to the proximity of their real parts to the zero eigenvalue (which is doubly degenerate —for the identity and the swapping operator). We consider only one of the two complex eigenvalues forming a conjugate pair. It is evident that all the gaps are closing with growing $N$. Eigenvalue crossings are also visible, while for moderate pump strengths there is an eigenvalue with an inverted-parabolic shape, flattening to zero as $N$ increases. This eigenvalue has zero imaginary part, suggesting that it is itself responsible for a first-order DPT in the thermodynamic limit mediated by the presence of bistability Minganti et al. (2018).

In addition, there are two eigenvalues with a nonzero imaginary part at low and high pump strengths. Their imaginary parts converge very quickly for relatively small $N$ (at $\tilde{F}/\gamma=1.65$ we have a maximum of $\sim 13$ bosons for $N=4$). These are the eigenvalues responsible for the limit cycles as $N\to\infty$. To illustrate this point, we show the frequencies of the semiclassical limit cycles (in black squares) in the plot $N=4$. For low pump values, the scaled frequency $\text{Im}[\lambda]/\gamma$ has already reached its $N\to\infty$ value for $N=4$. As we have already anticipated, two limit cycles with different periods coexist in the semiclassical bistable region ($1.02\lesssim\tilde{F}/\gamma\lesssim 1.1$). Either one or the other is reached in the long-time limit, depending on the amplitude of the initial coherent state.

The closure of the first gap for all the pump values we have considered suggests a critical slowing down transcending the bistability region, in contrast with local dissipation Note (2). To show this we employ the time-delayed second-order correlation function Fink et al. (2018) (symmetric with respect to the two modes $i=1,2$) $g^{(2)}(\tau)=\braket{\hat{a}_{i}^{\dagger}(0)\hat{a}_{i}^{\dagger}(\tau)\hat{a}_{i}(\tau)\hat{a}_{i}(0)}_{\rm ss}/\displaystyle\braket{\hat{a}_{i}^{\dagger}\hat{a}_{i}}_{\rm ss}^{2}=\operatorname{Tr}[\hat{a}_{i}^{\dagger}(0)\hat{a}_{i}(0)e^{\tau\mathcal{L}}(\hat{\rho}^{\prime})]/\braket{\hat{a}_{i}^{\dagger}\hat{a}_{i}}_{\rm ss}$, where $\hat{\rho}^{\prime}=\hat{a}_{i}\hat{\rho}_{\text{ss}}\hat{a}_{i}^{\dagger}/\operatorname{Tr}[\hat{a}_{i}^{\dagger}\hat{a}_{i}\hat{\rho}_{\text{ss}}]$ is the state of the system after an initial detection of the mode $i$ and is evolved during a time $\tau$ before making a second measurement. The subscript ${\rm ss}$ refers to the steady state. After a sufficiently long time, if the steady state is unique, $g^{(2)}(\tau)$ will relax to the value of $1$, due to its normalization. In Fig. 1(c) we plot $g^{(2)}(\tau)$ for three pump values. One can observe that the correlations take a long time to decay for $\tilde{F}/\gamma=0.5$, far below the bistability region. Furthermore, the function does not always converge to $1$. This is a clear manifestation of the infinite number of steady states: The initial detection changes the value of the conserved quantity $\langle\hat{Z}_{2}\rangle$, so the system must relax to a different steady state.

Figure 2 clearly depicts the emergence of the periodic oscillations as $N\to\infty$. In frames (b) and (c), we can observe frequency bistability, starting from two different initial conditions. Here we are in the semiclassical bistability region with finite $N$, where the oscillations are damped. Moreover, we show in Ref. Note (2) that single trajectories Fourier transform as well to a neat frequency peak matching that of the semiclassical response (apart from a small deviation due to stochastic jumps).

In the bistability region, a single quantum trajectory evidences switching between the dim and bright metastable states Carmichael (2015). This is accompanied by a frequency switch for the oscillations between the two (frequency) branches drawn in Fig. 1(c) for $N=4$. However, for large but finite $N$, the inverse of the switching rate is larger than the lifetime of the oscillations, which means that on average one obtains an oscillation frequency corresponding to either the lower or the upper branches of Fig. 1(c), depending on the initial condition [see Figs. 2(b) and 2(c)].

In Fig. 2(c), one can notice a mismatch between the frequency peaks in the semiclassical and the quantum regimes for $N=25$. There are two reasons for this: (i) the semiclassical evolution remains transient and the frequency has not yet attained its limit-cycle value (which depends nonlinearly on the mode population and therefore decreases in the long-time limit), (ii) the imaginary part of the Liouvillian eigenvalue responsible for the oscillations has not yet converged to its $N\to\infty$ value.

We now propose a mechanism accounting for the observed periodic oscillations. When $\tilde{U}=0$, bonding and antibonding modes are decoupled, and the antibonding mode evolves coherently giving trivial periodic oscillations in both the quantum and semiclassical regimes. When $\tilde{U}\neq 0$, the interaction term becomes very inefficient in coupling both modes when the antibonding-mode population is sufficiently small. The emergence of the limit cycles is then due to an effective decoupling between the two modes (see also Ref. Note (2) for further details).

IV Discussion

Nonlocal dissipation for open quantum systems without detailed balance has been considered in the context of noiseless subsystems or decoherence-free subspaces Zanardi and Rasetti (1997); *Zanardi98; Lidar et al. (1998), in cold atoms or many-body spin systems Iemini et al. (2018); James (1993); Olmos et al. (2014); Parmee and Cooper (2018), and in a recent proposal of a Floquet time crystal Wang et al. (2018). In general, that type of dissipation is either engineered or appears as a consequence of photon-mediated long range interactions. The specific form of Eq. (1) can be derived Note (2) from first principles if one considers that both modes in our system have the exact same coupling with a single bosonic bath. This approach is similar to Ref. Zanardi and Rasetti (1997); *Zanardi98.

In practice, the assumption of identical coupling between the two modes and the environment is idealistic. While the coupling strength between the two system modes and the reservoir modes (with frequencies in the vicinity of the driving frequency) could be equal, in general the relative phases will be different Note (2). This implies a mixture of bonding and antibonding dissipation and, in turn, the mixture of local and nonlocal dissipation following a change of basis. An important question, then, concerns the implications of a broken symmetry with respect to the swapping operation in the presence of a small yet non-negligible amount of dissipation for the antibonding mode. This perturbation has the consequence of lifting the degeneracy of the zero eigenvalue in the Liouvillian spectrum and also destroying the limit cycles. However, as we show in Ref. Note (2), periodic oscillations persist for long times, and the time-crystalline period is robust.

The BHD has already been studied in the laboratory Rodriguez et al. (2017). The type of dissipation in our model could be realized in a circuit QED experiment coupling two degenerate modes to a microwave resonator at a single spatial location, where the wave function overlap between the resonator mode and both system modes is identical. Alternatively, it could be realized with polaritons in two coupled semiconductor micropillars if the bonding and antibonding mode linewidths differ by two orders of magnitude or more. This could be achieved by etching the sample in such a way that the nonradiative losses Zambon et al. (2019) of the bonding mode make the dominant part of the total dissipation.

Acknowledgements.

We are thankful to the developers of the Quantum Toolbox in Python (QuTiP) Johansson et al. (2012, 2013), as we have used it for most of our calculations in the quantum regime. We thank A. Amo, E. Ginossar, and M. Stern for helpful discussions, and J. Keeling for the careful reading of the manuscript. C.L. gratefully acknowledges the financial support of CONICYT through Becas Chile 2017, Contract No. 72180352. M.H.S. gratefully acknowledges financial support from QuantERA InterPol and EPSRC (Grant No. EP/R04399X/1 and No. EP/K003623/2).