Steady states of a driven dissipative dipolar XXZ chain

TL;DR

This paper investigates the steady states of a driven dissipative dipolar XXZ spin chain, revealing a complex phase diagram with uniform, non-uniform, and bistable phases through mean-field and exact diagonalisation methods.

Contribution

It introduces a comprehensive analysis combining mean-field and exact diagonalisation to identify diverse phases and bistabilities in a driven dissipative dipolar XXZ chain.

Findings

Rich phase diagram with multiple phases

Fluctuations reveal non-uniform phases

Bistability signatures are general features

Abstract

We study theoretically a driven dissipative one-dimensional XXZ spin$-1/2$ chain with dipole coupling and a tunable strength of the Ising and XY interaction. Within a mean-field approximation, we find a rich phase diagram with uniform, spin density wave, antiferromagnetic and oscillatory phases, as well as regions of phase bistability. We study the phase diagram of small quantum systems using exact diagonalisation, and compare the results to the mean-field theory. We find that while expectation values only capture the uniform phases of the mean-field theory, fluctuations about these expectation values give signatures of spatially non-uniform phases and bistabilities. We find these signatures for all ratios of the Ising to XY interaction, showing that they appear to be general features of spin$-1/2$ systems

Equations

$$\begin{array}{rl}\frac{d\hat{\rho }(t)}{dt}=& -\text{i}\left[\hat{H},\hat{\rho }(t)\right]\\ & +\frac{\Gamma }{2}\sum _i^N\left[2{\hat{\sigma }}_i^{-}\hat{\rho }(t){\hat{\sigma }}_i^{+}-\left\{{\hat{\sigma }}_i^{+}{\hat{\sigma }}_i^{-},\hat{\rho }(t)\right\}\right],\end{array}$$

$$\begin{array}{rl}& \hat{H}=\frac{\Omega }{2}\sum _i^N{\hat{\sigma }}_i^x-\frac{\Delta }{2}\sum _i^N{\hat{\sigma }}_i^z\\ & +\sum _{i\ne j}^N\frac{J}{4R_{ij}^3}\left[\cos \alpha {\hat{\sigma }}_i^z{\hat{\sigma }}_j^z+\sin \alpha \left({\hat{\sigma }}_i^y{\hat{\sigma }}_j^y+{\hat{\sigma }}_i^x{\hat{\sigma }}_j^x\right)\right],\end{array}$$

$$\begin{array}{rl}\frac{d{S}_j^x}{dt}=& -\Gamma S_j^x-\Delta S_j^y\\ & +2\sin \alpha \sum _{i(\ne j)}^N\frac{J}{R_{ij}^3}{S}_j^z{S}_i^y-2\cos \alpha \sum _{i(\ne j)}^N\frac{J}{R_{ij}^3}{S}_i^z{S}_j^y,\\ \frac{d{S}_j^y}{dt}=& -\Gamma S_j^y-\Omega S_j^z+\Delta S_j^x\\ & -2\sin \alpha \sum _{i(\ne j)}^N\frac{J}{R_{ij}^3}{S}_j^z{S}_i^x+2\cos \alpha \sum _{i(\ne j)}^N\frac{J}{R_{ij}^3}{S}_i^z{S}_j^x,\\ \frac{d{S}_j^z}{dt}=& -\Gamma \left(S_j^z+\frac{1}{2}\right)+\Omega S_j^y\\ & -2\sin \alpha \sum _{i(\ne j)}^N\frac{J}{R_{ij}^3}(S_i^y{S}_j^x-S_j^y{S}_i^x),\end{array}$$

$$\begin{array}{r}\text{IoD}=\frac{{\sum }_{i,j}^N\left(⟨{\hat{S}}_i^z{\hat{S}}_j^z⟩-⟨{\hat{S}}_i^z⟩⟨{\hat{S}}_j^z⟩\right)}{{\sum }_i^N\left(⟨{\hat{S}}_i^z⟩+1/2\right)},\end{array}$$

$$\begin{array}{rl}\hat{\rho }(t)& =\exp (\mathcal{L}t)\hat{\rho }(0)\\ & ={\hat{\rho }}_{ss}+\sum _{i=1}^{4^N-1}\text{tr}{\hat{L}}_i\hat{\rho }(0){\hat{R}}_i{e}^{{\lambda }_{i}t},\end{array}$$

$$\begin{array}{rl}& ⟨\hat{A}(t+\tau )\hat{B}(t)⟩=\text{Tr}\{\hat{A}{e}^{(\tau +t-t)\mathcal{L}}\hat{B}\hat{\rho }(t)\}\\ & =\text{Tr}\{\hat{B}\hat{\rho }(t)\}\text{Tr}\{\hat{A}{\hat{\rho }}_{ss}\}+\sum _{i=1}^{4^N-1}{e}^{{\lambda }_i\tau }\text{Tr}\{{\hat{L}}_i\hat{B}\hat{\rho }(t)\}\text{Tr}\{\hat{A}{\hat{R}}_i\},\end{array}$$

$$\begin{array}{rl}& \underset{t\to \infty }{\lim }⟨\hat{A}(t+\tau )\hat{B}(t)⟩\\ & =\text{Tr}\{\hat{A}{\hat{\rho }}_{ss}\}\text{Tr}\{\hat{B}{\hat{\rho }}_{ss}\}+\sum _{i=1}^{4^N-1}\text{Tr}\{{\hat{L}}_i\hat{B}{\hat{\rho }}_{ss}\}\text{Tr}\{\hat{A}{\hat{R}}_i\}{e}^{{\lambda }_i\tau }\\ & \approx \text{Tr}\{\hat{A}{\hat{\rho }}_{ss}\}\text{Tr}\{\hat{B}{\hat{\rho }}_{ss}\}+\text{Tr}\{{\hat{L}}_1\hat{B}{\hat{\rho }}_{ss}\}\text{Tr}\{\hat{A}{\hat{R}}_1\}{e}^{{\lambda }_1\tau },\end{array}$$

Full text

I Introduction

The study of quantum systems driven far from equilibrium has attracted much interest over the last few years. Although not as well understood as their equilibrium counterparts, non-equilibrium phenomena are in fact rather prevalent, as any experiment has some form of interaction with an environment, which will induce dissipative processes. Whereas these processes could be viewed as a nuisance, recent studies have shown that the interplay between an external drive and dissipation can produce exotic non-equilibrium phases such as spin density waves (SDW), antiferromagnetic phases (AFM), persistent long-time oscillations (OSC) and phase bistabilities within spin$-1/2$ systems Chan et al. (2015); Lee et al. (2013, 2011); Wilson et al. (2016); Qian et al. (2015); Parmee and Cooper (2018) and higher spin systems Qian et al. (2012). Therefore, understanding the long-time steady state phases that can occur in open quantum systems is an intense area of current research.

In order to capture the steady state phases in a macroscopically large system, it is common to employ a mean-field approximation, where spatial correlations between sites are ignored. While the use of the mean-field approximation is well understood in equilibrium phenomena, an open question that still remains is the validity of the mean-field approximation for dissipative systems. For thermal gases in three-dimensions, one expects the mean-field approximation to become valid as quantum fluctuations become negligible Carr et al. (2013); Šibalić et al. (2016). However, for cold systems and/or systems in reduced dimensions, where quantum fluctuations are more important, this is not necessarily the case. There have already been many studies into the true phases and transitions for the driven dissipative spin systems. These have involved exact diagonalisation and quantum Monte Carlo wavefunction approaches for small system sizes Olmos et al. (2014); Rota et al. (2018); Lee et al. (2012); Ates et al. (2012), or incorporating the use of Keldysh methods Maghrebi and Gorshkov (2016), cluster mean-field Jin et al. (2016, 2018); Owen et al. (2018); Biella et al. (2018), Gutzwiller Monte Carlo approaches Huybrechts and Wouters (2019), corner-space renormalisation Rota et al. (2017) or variational approaches Weimer (2015) and Matrix Product and tensor network methods Cui et al. (2015); Mascarenhas et al. (2015); Mendoza-Arenas et al. (2016); Joshi et al. (2013); Kshetrimayum et al. (2017); Höning et al. (2013) to study larger systems. These studies have shown that first order transitions in the mean-field approximation can become second order when quantum fluctuations are included, and that bistabilities can be lost. They have also shown that the emergence of certain phases, such as antiferromagnetic phases or long-time oscillations may not occur in low dimensional spin-$1/2$ systems.

Despite the disagreements between mean-field and exact numerics of quantum systems, the mean-field approximation can still serve as an indicator of features that emerge in the full quantum system. For example, it has been shown that regions of mean-field bistability for the Ising model correspond to long spatial correlations in the full quantum model Hu et al. (2013). Also, while bistability has not been observed in finite-sized quantum spin systems, the bistable nature of the mean-field solutions is evident in quantum trajectories of the system Wilson et al. (2016); Lee et al. (2012); Olmos et al. (2014); Ates et al. (2012) and also results in a decrease in the spectral gap of the Liouvillian Wilson et al. (2016). Therefore, it is interesting to ask how general these features are when comparing quantum results to those within the mean-field approximation, and if there are other consequences of the mean-field results for the full quantum dynamics.

In this paper, we study a driven-dissipative XXZ model with a tunable XY and Ising interaction and dipole coupling as a function of detuning and external drive strength. The tunable Ising and XY interaction connects the Ising model in Ref. Lee et al. (2011) and the XY model in Ref. Wilson et al. (2016). By calculating the nonequilibrium phase diagram at mean-field level, we find how the phases evolve between these two limits, finding the emergence of spin density waves, antiferromagnetism, temporal oscillations and bistabilities. We then analyse small quantum systems and carry out an in-depth comparison to our mean-field phase diagram. While we do not find the same phase diagram for the full quantum system, we do find strong signatures of the mean-field in the spin fluctuations which relate to bistabilities and spatial phases that arise in the mean-field results. The tunable Ising to XY interaction allows us to track these signatures between the XY and Ising limit, and show that the quantum results are general features.

The rest of this paper is organised as follows. In Sec. II, we describe the model, then in Sec. III, we derive the mean-field phase diagram. In Sec. IV we look at the full quantum model and discuss the results in Sec. V before drawing conclusions in Sec. VI.

II Model

The system consists of a large number, $N$, of atoms or polar molecules in a one-dimensional (1D) array. Two internal dipolar energy levels, $\ket{g}$ and $\ket{e}$, are isolated and the transition between them driven by an external drive detuned from resonance. The system can then be modelled as a spin$-1/2$ system, interacting via dipole-dipole interactions with a Rabi drive and detuning. The two-level transition also has a finite lifetime, comparable to the time scales of interaction, which results in decay from the excited state. Under the Born and Markov approximations, the dynamics of the system are described by the following master equation ($\hbar=1$)

[TABLE]

where the spin operators are defined as $\hat{\sigma}^{z}_{i}=\ket{e_{i}}\bra{e_{i}}-\ket{g_{i}}\bra{g_{i}}$, $\hat{\sigma}^{-}_{i}=\ket{g_{i}}\bra{e_{i}}$ and $\hat{\sigma}^{+}_{i}=\ket{e_{i}}\bra{g_{i}}$, with $\ket{e_{i}}$ and $\ket{g_{i}}$ being the excited and ground states of the two-level transition respectively on site $i$. The decay constant, $\Gamma$, is controlled by optical pumping Lee et al. (2013), where the two level system is off-resonantly dressed by a third energy level. The effective decay constant of the two-level transition is then a combination of the decay rate of the original two-level system and the decay rate of the third energy level. This allows us to tune the decay independently from the interaction between the two-level systems. The Hamiltonian is given by

[TABLE]

where $R_{ij}=|r_{i}-r_{j}|$, with $r_{i}$ being the position vector of a two-level system at site $i$, and $\hat{\sigma}^{x}_{i}=\ket{e_{i}}\bra{g_{i}}+\ket{g_{i}}\bra{e_{i}}$ and $\hat{\sigma}^{y}_{i}=\text{i}\ket{g_{i}}\bra{e_{i}}-\text{i}\ket{e_{i}}\bra{g_{i}}$. The detuning is given by $\Delta=\omega_{d}-\omega_{eg}$, with $\omega_{d}$ being the drive frequency and $\omega_{eg}$ the two-level transition frequency, and the drive strength is given by the Rabi coupling $\Omega$. We note that the power law coupling is motivated by experiment but is not crucial for the results presented here, i.e. a nearest-neighbour coupling would also work. The dipole-dipole interaction is given by $J=J_{0}a^{3}(1-3\sin^{2}\Theta)$ where $a$ is the lattice spacing (see Fig. 1) and the parameter $J_{0}$ is given by $J_{0}=|\textbf{d}|^{2}/4\pi\epsilon_{0}a^{3}$ .

The angle $\Theta$ is the orientation of the dipoles, and can be tuned by application of a d.c. electric field. The parameter $\alpha$ relates the relative strength of the Ising and XY dipole interactions and can take values between $-\pi$ to $\pi$. However, we only focus on the range $0\leq\alpha<\pi$ as values below zero simply corresponds to a change in the sign of $J$. Tuning the value of $\alpha$ depends on the choice of internal states and external fields, with $\cos\alpha$ being related to the difference in dipole moments of the groundstate and excited state, and $\sin\alpha$ being related to the transition dipole moment between the groundstate and excited state Gorshkov et al. (2011); Peter (2015).

For the remainder of the paper, we work with $\Theta=\pi/2$. Other values of $\Theta$ will result in a sign change and scaling of the interaction in 1D, but will not lead to significant changes in the types of phases that appear in our system, only in the size of the regions as a function of $\Delta$ and $\Omega$. Also, we will only study $0\leq\alpha<\pi$, as the values $-\pi\leq\alpha<0$ will result in the same phases as for $0\leq\alpha<\pi$ but with the phase diagram reflected about $\Delta=0$ due to the sign change in the interaction term. Therefore, it is sufficient to consider the range $0\leq\alpha<\pi$ and $\Theta=\pi/2$ to cover all phases that can occur in the 1D system. Finally, we choose a value of $J_{0}a^{3}/\Gamma=5$ for the rest of this paper. We should find similar results for other nearby values of $J_{0}a^{3}/\Gamma$, but expect that for $J_{0}a^{3}/\Gamma\ll 1$, we will only find spatially uniform phases in the system as the spins become effectively decoupled.

In order to find the long-time steady state of Eq. (1), we make a Gutzwiller mean-field approximation which results in taking $\hat{\rho}(t)=\otimes\hat{\rho}_{i}(t)$, where $\hat{\rho}_{i}(t)$ is a density matrix on site $i$, and we effectively ignore correlations between spins. Then, by taking the trace of Eq. (1) over all the sites except a given site $j$, we obtain the equations of motion as

[TABLE]

where $S^{\beta}_{j}=\frac{1}{2}\Tr\left[\hat{\sigma}^{\beta}_{j}\hat{\rho}(t)\right]$ are the spin expectation values. To find the steady state solutions, we solve the dynamics of the non-linear Eqs. (3) by evolving them into the long-time limit.

III Mean-Field Phase Diagram

We compute the mean-field phase diagram by finding the steady states of Eqs. (3). To do this, we first employ a bipartite sublattice ansatz to find the uniform and antiferromagnetic steady state solutions. Depending on the parameters, there can be up to three uniform solutions and three sets of antiferromagnetic solutions. To determine the final phases that exist in the system, we perform linear stability analysis of the resultant solutions to fluctuations with wave vectors $ka=\pi m/N$, where $N$ is the number of sites on the lattice and $m$ is an integer in the range $0\leq m\leq N$. In cases where the wave vector of instability is not equal to [math] or $\pi$, we expect spin density wave solutions to form and for the bipartite sublattice ansatz to fail. To confirm the sublattice ansatz results are correct, we solve the full dynamics of Eqs. (3) by evolving the equations in time until the long-time limit (up to $t\Gamma=200$) for a system size of $N=100$ with periodic boundary conditions. Simulating the full dynamics also allows us to find the resultant phases in regimes where the sublattice ansatz breaks down. For our dynamical simulations, as an initial condition, we use either $(S^{x},S^{y},S^{z})=(0,0,-1/2)$ or, if examining phase instability, the steady state that becomes unstable. In Fig. 2, we show a collection of phase diagrams as a function of $\Delta$ and $\Omega$ for select values of $\alpha$ in the range $0\leq\alpha<\pi$.

The phase diagram for $\alpha=0$ is similar to a nearest-neighbour Ising model studied in Ref. Lee et al. (2011) as expected. Similarly, the phase diagram for $\alpha=\pi/2$ is similar to a nearest-neighbour XY model studied in Ref. Wilson et al. (2016). However, our phase diagrams now show how these phases change as the $\alpha$ value is moved away from these two cases. For all, or almost all, $\alpha$ values, we can see some general features that occur. Specifically, we can classify four key phases that emerge in the system. Firstly, for all $\alpha$, there are the spatially uniform phases, which are shown by the white regions in the phase diagrams. At low Rabi drive, the uniform phase has a high spin magnitude and the spins close to the state with $S^{z}\approx-1/2$. We define this uniform phase as the $U_{1}$ phase. At high drive, the spin magnitude decreases and the spins lie close to the state with $S^{z}\approx 0$. We denote this as the $U_{2}$ phase. Both uniform phases occur for any spin half system even without interactions as they are just general solutions to the optical Bloch equations. For $\alpha=\pi/4$, we have a Heisenberg Hamiltonian which conserves total spin and results in only uniform phases. For all $\alpha$, the $U_{1}$ phase smoothly crosses over into the $U_{2}$ phase for most parameter ranges. However, for $\alpha\neq\pi/4$, when the drive and detuning are comparable to the interaction strength ($J/\Gamma=-5$), regimes of bistability between the $U_{1}$ and $U_{2}$ exist, which lead to sharp transitions between the two phases. Which phase the system ends up in within this region depends on the initial conditions. These regions of bistability are denoted by the dark blue regions in the phase diagram.

It is also the case that when the drive and detuning are comparable to the interactions, for all $\alpha\neq\pi/4$, the uniform phases can become unstable to fluctuations, breaking translational invariance and giving rise to non-trivial phases. In the red regions, the uniform phase becomes unstable to fluctuations with a wavevector of $ka=\pi$ and a stable set of antiferromagnetic solutions exist. This results in the emergence of a canted AFM solution, with the nature of the AFM phases depending on the $\alpha$ value. The AFM solutions have the largest deviation between the $S^{z}$ components when $0<\alpha<\pi/4$, and a large deviation in the $S^{x}$ and $S^{z}$ components when $3\pi/4<\alpha<\pi$. However, when the XY interaction dominates, the AFM solution has the strongest deviation in the $S^{y}$ components. As well as instabilities to $ka=\pi$, the uniform phases can become unstable to $ka<\pi$ which results in the emergence of a SDW phase. In the SDW phase, shown by the green regions, the spin orientation varies periodically through the lattice with a period set by the instability wavevector, $ka$. We find that there are no SDW instabilities for $\alpha=0$. Therefore, it seems that SDW instabilities are related to presence of the XY interaction. The final key phase to emerge are persistent long-time oscillations, denoted by OSC, where the effects of the drive dominate over the effects from dissipation. We find that all the oscillations emerge from the instability of the AFM phase, which undergoes a Hopf bifurcation, and inherit an AFM nature. The oscillations occur in the pink regions of the phase diagrams.

In several regions of the phase diagrams, multiple solutions to Eqs. (3) coexist, which can lead to bistabilities. In the yellow regions, a stable uniform solution and stable set of AFM solutions exist, which results in AFM/$U_{1}$ bistability. Similarly, in the light yellow regions, there is an OSC phase which is also bistable to the $U_{1}$ phase. We do find that there are cases where the oscillations become unstable and so only the uniform phase exists, which we have not marked in our phase diagram. In the light blue regions, both the $U_{1}$ phase and $U_{2}$ phase exist, but one becomes unstable to fluctuations. When simulating the dynamics in these regions, we find there is predominantly only one stable uniform solution and only small regions of $U_{2}$/SDW phase bistability. Similarly, in the orange regions, there exists both a uniform phase and an antiferromagnetic phase, but the uniform phase is unstable. We find when simulating the dynamics in this region that only the AFM phase exists and that the AFM phase is frustrated, with the amount of frustration depending on the initial conditions. Likewise, in the light orange regions, there is an OSC phase and an unstable uniform phase. We find only the OSC phase exists, and again that there are frustration effects.

IV Quantum Phase Diagram

Having now calculated the phase diagram at mean-field level, we now examine how it compares to the phase diagram of the full quantum system. To do this, we look at the steady state of small quantum systems with $N=8$ spins and periodic boundary conditions. Despite the small system size, we can see some distinct features emerge for the quantum system that reflect the mean-field phases.

We first examine the spin expectation values on each site for a direct comparison to the mean-field results. We find that for the quantum system, the expectation values on each site are uniform, with no spatial variation, for all values of detuning and drive. Even for small systems, we would expect some non-uniform phases to emerge within the mean-field approximation, as the wave vectors that cause instability of the uniform phases, such as $ka=\pi$, are permitted. However, we find only uniform phases for the small quantum system, which we would expect to persist until large system sizes, due to the translational invariance of the system. If we compare the quantum expectation values to the mean-field expectation values from our phase diagrams in the regions where stable uniform phases exist, we find there is good agreement when the magnitude of the Rabi drive and detuning are large, for all $\alpha$. The difference between the expectation values becomes larger in regions where there is $U_{1}$/$U_{2}$ bistability. This is because we find no bistability in small quantum systems, but a smooth crossover between the $U_{1}$ and $U_{2}$ phases. Therefore the quantum expectation values will eventually differ from either choice of the $U_{1}$ or $U_{2}$ mean-field solution. For $\alpha=\pi/4$, we find exact agreement between the quantum and mean-field expectation values for all detuning and Rabi drive due to the Heisenberg symmetry.

Although the expectation values only show a single uniform phase, the connected correlators between sites give insight into spatial structure of the fluctuations about the expectation values and possible emergence of non-uniform phases in large systems. In Fig. 3 (a), we plot the connected correlator, $\langle\hat{S}^{z}_{i}\hat{S}^{z}_{j}\rangle_{c}=\langle\hat{S}^{z}_{i}\hat{S}^{z}_{j}\rangle-\langle\hat{S}^{z}_{i}\rangle\langle\hat{S}^{z}_{j}\rangle$, between a site $i$ and its nearest-neighbour for $\alpha=0$. We find that the connected correlator changes sign between nearest-neighbour sites in the region where where AFM solutions exist in the mean-field, but maintains the same sign when in the uniform region. Figure 3 (b) and (c) show examples of how the connected correlator varies across the lattice sites for a choice of $\Delta$ and $\Omega$ in the AFM region and in the uniform region, shown by the red and blue circles respectively. We can see how both connected correlators lose long-range order quickly, but maintain an alternating sign in the AFM region, while being persistently positive in the uniform region. We also plot the $\langle\hat{S}^{y}_{i}\hat{S}^{y}_{j}\rangle_{c}$ connected correlator for $\alpha=\pi/2$ in Fig. 4, finding again a loss of long-range order, but sign changes between nearest-neighbour sites in the region where AFM solutions exist in the mean-field phase diagram. The $\alpha=0$ connected correlator results are similar to those in Ref. Lee et al. (2011) and the $\alpha=\pi/2$ results similar to those in Ref. Wilson et al. (2016). For both Fig. 3 and Fig. 4, the choice of connected correlator is based on the spin component with the strongest deviation from the spatially uniform phase in the mean-field analysis.

We find the strongest signature of the mean-field occurs when looking at phase bistability. In regions of bistability, many studies Wilson et al. (2016); Olmos et al. (2014); Lee et al. (2012) have found examples of bimodality in the expectation value distributions when using Quantum Monte Carlo wavefunction methods. This bimodality arises when there are two stable mean-field solutions where, in the mean-field limit, each steady state behaves as if the minimum of a double well potential. Classically, the system will sit in either steady state indefinitely. However, in the full quantum system, quantum fluctuations induce transitions between these steady states, giving a single expectation value, but leaving a double well structure in the expectation value distributions. For the phase diagram in Fig. 2, bistability brings the greatest change in $S^{z}$ values between the two mean-field phases, and so, for the quantum system, we expect to see bimodality in the total number of excitations, $\sum_{i}^{N}(\hat{S}^{z}_{i}+\hat{\mathbb{1}}/2)$. As first studied in Ref. Wilson et al. (2016), a bimodal distribution can be observed by studying the Index of Dispersion (IoD), which is given by

[TABLE]

and is a measure of a distribution’s variance normalised by its mean. In the limit of zero Rabi drive, when $\langle\hat{S}_{i}^{z}\rangle=-1/2$, we have $\text{IoD}=1$, whereas in the limit of high Rabi drive, when $\langle\hat{S}_{i}^{z}\rangle\approx 0$, we have $\text{IoD}=1/2$. In between these two limits, the IoD will either decrease or, when there is a bimodality in the distribution, increase above unity. In Fig. 5, we plot the IoD as a function of $\Delta$ and $\Omega$.

We can see from Fig. 5 that the agreement between the fluctuation region and the bistability region is good for all $\alpha$ values, suggesting that the large peak in the IoD is indeed due to the double well structure of the mean-field solution appearing in the quantum steady state. The results for $\alpha=\pi/2$ are similar to those in Wilson et al. (2016) for the nearest-neighbour XY model, but the results for other $\alpha$ are new.

Another signature of bistability should also be evident in spectral gap of the Liouvillian. The master equation, Eq. (1), can be written as $d\hat{\rho}(t)/dt=\mathcal{L}\hat{\rho}(t)$, where $\mathcal{L}$ is the Liouvillian superoperator that determines the evolution of the system. By expanding the density matrix in eigenvectors of the Liouvillian superoperator, we can write a generic quantum state as

[TABLE]

where $\hat{L}_{i}$ and $\hat{R}_{i}$ are the left and right eigenvectors of $\mathcal{L}$ respectively, with complex eigenvalues, $\lambda_{i}$. The eigenvalues of the Liouvillian come in complex-conjugate pairs due to the Hermitian nature of the density matrix, with a negative real part corresponding to the decay of the eigenvalues and the imaginary part corresponding to coherences. We order the eigenvalues, $\lambda_{i}$, by the increasing magnitude of their real part. There is always at least one zero eigenvalue of the Liouvillian matrix ($\lambda_{0}=0$) corresponding to the steady state of the system, which we have extracted from the sum in Eq. (5) and denoted the eigenvector as $\hat{\rho}_{ss}\equiv\hat{R}_{0}$. Note that $\hat{\rho}_{ss}$ has unit trace and the remaining eigenvectors, $\hat{R}_{i}$, are traceless.

In systems with bistability, the spectral gap (i.e. the eigenvalue with smallest real part) should close so the system has two steady states. However, for small quantum systems, it is known that the gap will remain finite Schirmer and Wang (2010). Despite this, bistability may occur in very large quantum systems, and if this is to be the case, then we expect the gap to decrease in the region of bistability even for small systems. In Fig. 6, we compute Navarrete-Benlloch (2015) and plot the real part of the spectral gap for the system with $N=6$ spins.

We see there is indeed a reduction in the gap size in the bistability region compared to elsewhere in the phase diagram. We also study the real part of the spectral gap in smaller system sizes ($N=3,4,5$) and find that the real part of the spectral gap decreases with increasing $N$. This may indicate a closing of the gap in larger system sizes, although the gap could also saturate.

As well as examining the real part of the spectral gap, we also analyse the imaginary part. We expect that in regions where there are persistent oscillations, that the real part of the gap should go to zero and the imaginary part should have a non-zero value. However, we have found little sign this was the case and have not shown the data here. This behaviour may be due to the fact that the oscillations come from instabilities of the AFM phases. Our analysis of the real part of the spectral gap showed little decrease in the real part of the spectral gap in regions where there are mean-field bistabilities between the uniform and the AFM phases. Therefore, it appears to be the case that non-uniform phases give weak spectral signatures, which may explain the absence of any spectral signature for oscillations. It would be interesting to examine this further.

To measure the spectral gap in experiment, one could look at two-time correlators, whose decay depends specifically on the spectral gap in the long-time limit. For two operators, $\hat{A}$ and $\hat{B}$, the two-time correlator is given by Breuer and Petruccione (2007)

[TABLE]

where we have inserted the density matrix expansion given in Eq. (5). If we allow $t\rightarrow\infty$ so we study the two-time correlator in the steady state, then Eq. (6) simplifies to

[TABLE]

where in the last line we have assumed that $\tau$ is large and that the spectral gap, $\lambda_{1}$, is well separated from the rest of the Liouvillian spectrum. This shows how the two-time correlator will decay exponentially with a decay time set by the spectral gap. Therefore, if the spectral gap does decrease in a given parameter regime, we should find that any two-time correlator will have a much longer temporal decay than in regions where the spectral gap is large. In Fig. 7, we plot the difference in real part of the first and second eigenvalues of the Liouvillian to see if the spectral gap is well separated from the spectral bulk. We find for a dominant Ising interaction, there is a separation inside the bistable regime, so measuring the gap from the connected correlator decay should be possible. For a dominant XY interaction, the gap from the bulk is quite small and therefore measuring the gap from the two-time correlator will be harder.

V Discussion

We have studied the mean-field phase diagram of an driven-dissipative XXZ model with a tunable XY and Ising interaction, finding that the interplay between drive and dissipation leads to four key types of non-equilibrium phases. Specifically, we find uniform phases, AFM phases, SDW phases and OSC phases, as well as bistabilities between these phases, and have studied how these phases change with the tuning of the Ising to XY interaction. Such a system could be readily studied with Rydberg atoms Glaetzle et al. (2015), trapped ions Britton et al. (2012) or polar molecules Gorshkov et al. (2011); Wall et al. (2015), where dissipation can be controllably induced with optical pumping Lee et al. (2013).

Our in-depth study of the full quantum system for a small number of spins shows that the expectation values of the small quantum system agree well with the mean-field uniform phases at strong Rabi drive and detuning, but do not agree as well when the drive and detuning are comparable to the interaction strength. The biggest difference between the mean-field and quantum results is the small quantum system does not exhibit any bistabilities, OSC, AFM or SDW phases. However, analysis of connected correlators shows how fluctuations about these expectation values give signatures of the non-spatially uniform mean-field phases and indicate the possibility of agreement between the mean-field and quantum results for large enough systems. The strongest signatures of the mean-field can be found when looking at mean-field bistability, where we have shown there is good agreement between the regions of mean-field bistability and where the IoD of excitations increases above unity. Finally, we have shown a good agreement between the bistability region and a decrease in the real part of the spectral gap in the Liouvillian, which should be observable by a increase in the decay times of two-time correlators.

As mentioned in earlier sections, many of our results here agree with similar studies with nearest-neighbour interactions Lee et al. (2011); Wilson et al. (2016), both at the mean-field and quantum level, indicating that the $1/r^{3}$ power-law nature of the interactions have little effect on the resultant phases. Given this and the fact that our system is 1D, modelling the system with DMRG Schollwöck (2005) to achieve larger system sizes should be possible. It would be interesting to see if the quantum fluctuations grow stronger in larger systems and lead to the breaking of spatial uniformity we have seen in our analysis, resulting in a phase diagram closer to our mean-field results. It would also be interesting to study how the IoD peak grows with system size. When studying smaller system sizes, we find the maximum value of the IoD for each $\alpha$ value grows with system size. If the mean-field is expected to become increasingly correct with increased dimensionality and system size, then the IoD should to drop to zero as the connected correlators in Eq. (4) become zero. However, we do not expect the mean-field to become valid for a 1D system. Looking at how the IoD behaves in higher dimensionalities would also be interesting. Finally, carrying out finite-size scaling of the Liouvillian gap to larger system sizes than studied here may indicate if the gap will eventually close in large enough quantum systems.

VI Conclusions

We have studied the mean-field phase diagram of a driven-dissipative XXZ model with a tunable Ising to XY interaction. We find the emergence of four key types of phase: uniform phases, spin density wave phases, antiferromagnetism and oscillatory phases as well as phase bistabilities. We characterise how the nature of these phases change with the relative strength of the Ising to XY ratio. We then analyse the phases of the corresponding quantum system for small numbers of spins and find that the mean-field results correspond to key features in the quantum phase diagram. We find that peaks in the index of dispersion and decreases in the real part of the Liouvillian gap appear in regimes where the mean-field theory shows bistability. We also find signatures of spatially varying mean-field phases in the connected correlators computed for the full quantum problem. These features persist with changes in the relative strength of the Ising and XY interaction, indicating that they are generic features of spin$-1/2$ systems.

Statement of compliance with EPSRC policy framework on research data: All data accompanying this publication are directly available within the publication.

VII Acknowledgements

This work was supported by EPSRC Grant Nos. EP/K030094/1 and EP/P009565/1 and by the Simons Foundation.