Excitations and correlations in the driven-dissipative Bose-Hubbard model

TL;DR

This paper investigates the driven-dissipative Bose-Hubbard model using a field-theoretic approach, revealing complex phase behavior, excitations, and correlations influenced by driving and losses.

Contribution

It introduces a detailed analysis of the phase diagram and excitation regimes of the driven-dissipative Bose-Hubbard model using the Schwinger-Keldysh formalism.

Findings

Identification of bistability and phase regimes

Discovery of dynamically unstable, gapped, and gapless phases

Correlation length increases in the gapless regime

Abstract

Using a field-theoretic approach within the Schwinger-Keldysh formalism, we study a Bose-Hubbard model in the presence of a driving field and dissipation due to one-body losses. We recover the bistability diagram from the Gross-Pitaevski equation and analyze the different phases with respect to their elementary excitations and correlations. We find the low-density solution to be subdivided into a dynamically instable, a gapped, and a gapless regime. The correlations decay exponentially, but a substantial increase of correlation length marks the regime of gapless excitations.

Equations

$$H_{\mathrm{BH}}=-J\sum _{⟨i,j⟩}{b}_i^{†}{b}_j+\sum _i\left(\frac{U}{2}\sum _i{b}_i^{†}{b}_i^{†}{b}_i{b}_i+{\omega }_{\mathrm{c}}{b}_i^{†}{b}_i\right),$$

$$H_{\mathrm{drive}}(t)=\sum _i\left(j{b}_i^{†}{\mathrm{e}}^{i{\omega }_{\mathrm{p}}t}+j^{\ast }{b}_i{\mathrm{e}}^{-i{\omega }_{\mathrm{p}}t}\right).$$

$$H=$$

$$i{\partial }_t\rho =[H,\rho ]+\frac{i\gamma }{2}\sum _i\left(2b_i\rho b_i^{†}-b_i^{†}{b}_i\rho -\rho b_i^{†}{b}_i\right),$$

$$S$$

$$\displaystyle-\frac{i\gamma}{2}\left(2b_{i+}(t)b_{i-}^{*}(t)-b_{i+}^{*}(t)b_{i+}(t)-b_{i-}^{*}(t)b_{i-}(t)\right)\Big{]}+J\sum_{\langle i,j\rangle}\left(b_{i+}^{*}(t)b_{j+}(t)-b_{i-}^{*}(t)b_{j-}(t)\right)\Big{\}}.$$

$$S$$

$$\displaystyle+J\sum_{\langle i,j\rangle}\left(b_{i,{\rm c}}^{*}b_{j,{\rm q}}+b_{i,{\rm q}}^{*}b_{j,{\rm c}}\right)\Big{\}}.$$

$$\frac{\delta S}{\delta b_{i,\mathrm{c}}^{\ast }}=0\mathrm{and}\frac{\delta S}{\delta b_{i,\mathrm{q}}^{\ast }}=0,$$

$$b_{i,\mathrm{q}}=0,\mathrm{and}{b}_{i,\mathrm{c}}=\sqrt{2}b,$$

$$b\left[\mu -U|b{|}^2+zJ-\frac{i\gamma }{2}\right]-j=0.$$

$$n\left[(\mu -nU+zJ{)}^2+\frac{{\gamma }^2}{4}\right]=|j{|}^2,$$

$$b_{i,\mathrm{q}}(t)={\delta }_{i,\mathrm{q}}(t),\mathrm{and}{b}_{i,\mathrm{c}}(t)=\sqrt{2}b+{\delta }_{i,\mathrm{c}}(t).$$

$$S^{(2)}=$$

$$\displaystyle-Ub^{2}\delta_{i,{\rm q}}^{*}\delta_{i,{\rm c}}^{*}+\frac{i\gamma}{2}|\delta_{i,{\rm q}}|^{2}\Big{]}+$$

$$\displaystyle+J\sum_{\langle i,j\rangle}\left(\delta_{i,{\rm c}}^{*}\delta_{j,{\rm q}}+\delta_{i,{\rm q}}^{*}\delta_{j,{\rm c}}\right)\Bigg{\}}+{\rm c.c.}$$

$$S^{(2)}=$$

$$\left(\begin{array}{cc}0& P^{\mathrm{A}}(\omega ,\mathbf{k})\\ P^{\mathrm{R}}(\omega ,\mathbf{k})& P^{\mathrm{K}}(\omega ,\mathbf{k})\end{array}\right)\cdot \left(\begin{array}{c}{\mathbf{\delta }}_{\mathrm{c},\mathbf{k}}(\omega )\\ {\mathbf{\delta }}_{\mathrm{q},\mathbf{k}}(\omega )\end{array}\right)$$

$$G^R(\omega ,\mathbf{k})=G^A(\omega ,\mathbf{k}{)}^{†}=$$

$$-i\left(\begin{array}{cc}⟨{\delta }_{\mathrm{c},\mathbf{k}}(\omega ){\delta }_{\mathrm{q},\mathbf{k}}(\omega {)}^{\ast }⟩& ⟨{\delta }_{\mathrm{c},\mathbf{k}}(\omega ){\delta }_{\mathrm{q},-\mathbf{k}}(\omega )⟩\\ ⟨{\delta }_{\mathrm{c},-\mathbf{k}}(\omega {)}^{\ast }{\delta }_{\mathrm{q},\mathbf{k}}(\omega {)}^{\ast }⟩& ⟨{\delta }_{\mathrm{c},-\mathbf{k}}(\omega {)}^{\ast }{\delta }_{\mathrm{q},-\mathbf{k}}(\omega )⟩\end{array}\right).$$

$$G^K(\omega ,\mathbf{k})=-G^A(\omega ,\mathbf{k})P^K(\omega ,\mathbf{k})G^R(\omega ,\mathbf{k}).$$

$$G^K(\omega ,\mathbf{k})=$$

$$-i\left(\begin{array}{cc}⟨{\delta }_{\mathrm{c},\mathbf{k}}(\omega ){\delta }_{\mathrm{c},\mathbf{k}}(\omega {)}^{\ast }⟩& ⟨{\delta }_{\mathrm{c},\mathbf{k}}(\omega ){\delta }_{\mathrm{c},-\mathbf{k}}(\omega )⟩\\ ⟨{\delta }_{\mathrm{c},-\mathbf{k}}(\omega {)}^{\ast }{\delta }_{\mathrm{c},\mathbf{k}}(\omega {)}^{\ast }⟩& ⟨{\delta }_{\mathrm{c},-\mathbf{k}}(\omega {)}^{\ast }{\delta }_{\mathrm{c},-\mathbf{k}}(\omega )⟩\end{array}\right).$$

$$P^R(\omega ,\mathbf{k})=-i\left(\begin{array}{cc}\omega +C_{\mathbf{k}}& -U{b}^2\\ -U{b}^{\ast 2}& -\omega +C_{\mathbf{k}}^{\ast }\end{array}\right),$$

$$P^K(\omega ,\mathbf{k})=i\gamma \left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right).$$

$${\omega }_{\pm }(\mathbf{k})=\pm \sqrt{(\mu +J_{\mathbf{k}}-2U|b{|}^2{)}^2-U^2|b{|}^4}-\frac{i\gamma }{2}.$$

$$n_{\mathbf{k}}\equiv ⟨{\delta }_{\mathbf{k}}{\delta }_{\mathbf{k}}^{\ast }⟩=\frac{i}{4\pi }\int G_{11}^K(\omega ,\mathbf{k})\mathrm{d}\omega -\frac{1}{2}.$$

$${\Delta }_{\mathbf{k}}\equiv ⟨{\delta }_{\mathbf{k}}{\delta }_{\mathbf{k}}⟩=\frac{i}{4\pi }\int G_{12}^K(\omega ,\mathbf{k})\mathrm{d}\omega ,$$

$${\Delta }_{\mathbf{k}}^{\ast }=⟨{\delta }_{\mathbf{k}}^{\ast }{\delta }_{\mathbf{k}}^{\ast }⟩=\frac{i}{4\pi }\int G_{21}^K(\omega ,\mathbf{k})\mathrm{d}\omega .$$

$${\tilde{n}}_0\equiv ⟨(b+{\delta }_0)(b^{\ast }+{\delta }_0^{\ast })⟩=|b{|}^2+n_0,$$

$${\tilde{\Delta }}_0\equiv ⟨(b+{\delta }_0)(b+{\delta }_0)⟩=|b{|}^2+{\Delta }_0,$$

$$⟨{\delta }_{\mathbf{r}}{\delta }_{\mathbf{r}+{\mathbf{r}}^{\prime }}^{\ast }⟩=\frac{1}{(2\pi {)}^2}\sum _{\mathbf{k}}{e}^{i\mathbf{k}\cdot (\mathbf{r}-{\mathbf{r}}^{\prime })}{n}_{\mathbf{k}}.$$

$$|⟨{\delta }_0{\delta }_x^{\ast }⟩|\propto \exp (-x/\xi ).$$

Full text

I Introduction

The Bose-Hubbard model is a paradigmatic model to describe a quantum phase transition, as site-to-site tunneling processes and on-site interaction processes compete energetically with each other fisher89 . In a groundbreaking experiment with cold atoms in an optical lattice, the Bose-Hubbard Hamiltonian was realized in a quantum simulator in 2002 greiner02 . Since then, there has been strong interest in exploiting light-matter interactions to implement also a photonic version of the Bose-Hubbard model hartmann06 ; greentree06 ; angelakis07 . The hopping term is then implemented as an array of coupled cavities, and the presence of either a real atom or an artificial atom in each cavity achieves the non-linearity of the model. A recent realization consists of a linear array of 72 cavities which are coupled to superconducting qubits fitzpatrick17 . Another realization of a Bose-Hubbard dimer, that is a system of two coupled cavities, has been described in Ref. rodriguez16 , using polaritons in a semiconductor to obtain the non-linear term. Due to the finite lifetime of the photons or polaritons, these systems do not conserve particle number. To avoid ending up in the vacuum, the losses have to be compensated by an appropriate driving term. A Bose-Hubbard model which is enriched by these two terms is usually refered to as the “driven-dissipative Bose-Hubbard model”, and has been studied in Refs. carusotto09 ; leboite13 ; leboite14 ; wilson16 ; cao16 ; casteels17 ; mamaev18 ; vicentini18 ; hannukainen18 ; biondi18 .

As the driving term explicitly breaks the U(1) symmetry, no quantum phase transition via spontaneous symmetry breaking occurs in the driven model. Nevertheless, the model is not featureless: In particular, it has been shown that the coupling between cavities induces a transition from a monostable to a bistable regime leboite13 ; leboite14 , although the variational principle used in Ref. weimer15 might suggest that the bistability is a mean-field artefact. On the other hand, renormalization group arguments for the XX model presented in Ref. maghrebi16 , as well as an interpretation of the full quantum dynamics in terms of non-equilibrium Langevin equations foss-feig17 , suggest that the driven-dissipative model falls into the same universality class as the equilibrium Ising model. In this light, the bistability is interpreted as the two minima of the Ising free energy, providing a non-equilibrium analog of the symmetry-broken phase.

In the present paper, we extend the mean-field treatment of the driven-dissipative Bose-Hubbard model, and calculate two-body correlators in order to distinguish between different parameter regimes. As a starting point, we employ a field-theoretic approach on the Keldysh time contour, which is suited to capture the dynamics of the corresponding master equation sieberer16 . By suppressing all fluctuations, this approach reproduces the results obtained from solving the dissipative Gross-Pitaevski equation, including also the bistability transition. We then admit for fluctuations up to quadratic order, and calculate the Green functions. With this, we gain new insight in the two different steady states. In particular, we calculate the momentum distribution, and we show that the low-density state has a peak structure different from the high-density phase. Studying the correlations in real space, we find an exponential decay in any regime, that is, no true long-range order is established in the system. However, an abrupt increase in the correlation length is seen when the hopping parameter is tuned. The division of the low-density phase into a short-ranged and a longer-ranged regime is also reflected by the elementary excitations of the system, as in the regime of increased correlation length the real part of the excitation spectrum vanishes at certain momenta.

The paper is organized as follows: Sec. II describes the model, Sec. III introduces the Schwinger-Keldysh method, and Sec. IV presents the results.

II Model

The Hamiltonian part of the system dynamics consists of the standard Bose-Hubbard model $H_{\rm BH}$ and a driving term $H_{\rm drive}(t)$:

[TABLE]

and

[TABLE]

Here, $J$ denotes the hopping strength, $U$ the interaction strength. The cavity frequency $\omega_{\rm c}$ acts like a chemical potential. The pump laser has an amplitude $j$ and a frequency $\omega_{\rm p}$. In a frame rotating with the pump frequency $\omega_{\rm p}$, the Hamiltonin becomes time-independent:

[TABLE]

Here, the on-site potential term is shifted by the pump frequency, $\mu\equiv\omega_{\rm p}-\omega_{\rm c}$.

Apart from the Hamiltonian dynamics, we will also take into account dissipative particle losses. The full system dynamics is then captured by a master equation

[TABLE]

with the dissipation rate $\gamma$.

III Schwinger-Keldysh action

In contrast to the time evolution of state vectors, the quantum-mechanical time evolution of an operator consists of a forward evolution and a backward evolution, $O(t)=U(t,0)^{\dagger}O(0)U(t,0)$. Therefore, any dynamics which goes beyond the evolution of pure state vectors requires treatment on a closed time contour. Accordingly, the Schwinger-Keldysh action is defined on two time branches, and for any quantity which enters the action we have to specify whether it lies on the forward or on the backward branch. Therefore, let us denote the complex-valued fields on the two branches by $b_{i+}(t)$ and $b_{i-}(t)$. The Schwinger-Keldysh action corresponding to Eq. (4) is written as

[TABLE]

It is convenient to rotate the fields to so-called classical fields $b_{i,{\rm c}}(t)=\frac{1}{\sqrt{2}}\left[b_{i+}(t)+b_{i-}(t)\right]$, and quantum fields $b_{i,{\rm q}}(t)=\frac{1}{\sqrt{2}}\left[b_{i+}(t)-b_{i-}(t)\right]$. In this basis, after integrating by parts a term containing $\partial_{t}b_{i,{\rm q}}$, the action reads

[TABLE]

For brevity, we have suppressed the time-dependence of the fields.

III.1 Mean-field solution

A mean-field solution to the equations of motion,

[TABLE]

is found by setting

[TABLE]

where $b$ has to fulfill a Gross-Pitaevski-like equation:

[TABLE]

The same equation has been formulated in Ref. leboite14 . Multiplied by its complex conjugate, it leads to a bistability criterion:

[TABLE]

where the mean photon number $n=|b|^{2}$ needs to be real. From this equation, a monostable phase is predicted when only one real solution exists, which is the case for small $J/\mu$. For larger $J/\mu$, Eq. (10) admits admits three real solutions. A comparison with the $P$-representation approach suggests that only two of them are stable leboite14 . Such comparison also suggests that the Gross-Pitaevski approximation underestimates the extension of the monostable regime leboite14 .

III.2 Fluctuations around mean-field

In the previous section, all fluctuations have been neglected. In the next step, we will take them into account up to second order. We write

[TABLE]

Via the choice of $b$ from Eq. (9), the first-order contribution to the action vanishes, $S^{(1)}=0$. Thus, to leading order we have an action which is quadratic in $\delta$:

[TABLE]

In general, the quadratic Schwinger-Keldysh action is given in terms of three Green functions. In Fourier space and after introducing Nambu spinors $\boldsymbol{\delta}_{{\rm c},{\bf k}}(\omega)=(\delta_{{\rm c},{\bf k}}(\omega),\delta_{{\rm c},{-\bf k}}(\omega)^{*})^{T}$, this relation reads

[TABLE]

Here, $P^{A}(\omega,{\bf k})=[G^{A}(\omega,{\bf k})]^{-1}$ and $P^{R}(\omega,{\bf k})=[G^{R}(\omega,{\bf k})]^{-1}$ are the inverse of the advanced and the retarded Green functions:

[TABLE]

From the retarded or advanced Green function, we obtain the excitation spectrum by solving for ${\rm det}[G^{R}(\omega,{\bf k})]^{-1}=0$. The function $P^{K}(\omega,{\bf k})$ is related to the Keldysh Green function $G^{K}$ in the following way:

[TABLE]

The Keldysh Green function contains information about correlation of the classical fields:

[TABLE]

Comparison of the general expression (III.2) with the Fourier-transformed action of Eq. (III.2) yields:

[TABLE]

where $C_{\bf k}=\mu+\frac{i\gamma}{2}-2U|b|^{2}+2J[\cos(k_{x})+\cos(k_{y})]$. The Keldysh component reads

[TABLE]

The excitation spectrum obtained from the retarded Green function matches with the spectrum found in Ref. leboite14 :

[TABLE]

From the diagonal elements of the Keldysh Green function, we obtain the momentum distribution:

[TABLE]

The off-diagonal elements are related to superfluid correlations:

[TABLE]

Note that both $n_{\bf k}$ and $\Delta_{\bf k}$ are defined as correlators between the fluctuating part $\delta$ of the fields. At ${\bf k}={\bf 0}$, we have a contribution from the mean-field part:

[TABLE]

and

[TABLE]

as $\langle\delta_{{\bf 0}}\rangle=0$.

IV Results

In the following section we will explicitly evaluate the action derived above, and extract relevant physical quantities. To ease our discussion, we will reduce the parameter space: Of the five energy scales ($\mu$, $U$, $zJ$, $\gamma$, and $j$), one scale can be taken as a unit of energy (here $\mu$), leaving four independent parameters. Since the truncation of the action to quadratic order can only be justified if the effect of interactions is small, we keep $U$ fixed at a small level, $U/\mu=0.5$. In Fig. 1, we have further set $\gamma$ equal to $j$, such that we can plot the two-dimensional parameter space. This plot shows the division into bistable and monostable phases according to Eq. (10), where the two solutions within the bistable regime can be distinguished by their density. Interestingly, the phase boundary of the bistable regime can also be obtained from the imaginary part of the excitation spectra, ${\rm Im}[\omega(\bf k)]$, which becomes gapless at the bistability transition. Along the lower branch of the boundary (from $B$ to $C$ in Fig. 1), the vanishing imaginary gap occurs to excitations above the low-density solution. Accordingly, the low-density solutions become instable here, and the system enters a monostable, high-density phase. The opposite is the case along the upper branch of the phase boundary (from $A$ to $C$), with a vanishing gap occuring in the imaginary part of the excitations above the high-density phase, and the system entering a monostable low-density phase.

The criterion of a vanishing gap in ${\rm Im}[\omega(\bf k)]$ also serves to establish a direct boundary between the two monostable phases. On such critical line (from $C$ to $D$ in Fig. 1), both types of excitations, high- and low-density excitations show not only gapless imaginary but also gapless real parts. Gapless real parts (without gapless imaginary parts) are also found within the bistable phase (red line).

IV.1 Low-density phase

We first study the elementary excitations obtained by plugging the low-density solution of Eq. (9) into Eq. (19). As illustrated in Fig. 2, the spectra vary significantly depending on the choice of parameters $J/\mu$ and $j/\mu$: If both cavity coupling $J$ and driving $j$ are weak, as illustrated in panel (a), the excitations are gapped. The gap can be removed by increasing either $J$ or $j$. As shown in panel (b) and (d), increasing $j$ leads to an extended region in $\bf k$-space in which the real part of the excitation frequency, ${\rm Re}(\omega)$, vanishes. In these regions, the imaginary part of the frequencies, gains a contribution additional to the otherwise constant value $-i\gamma/2$. This constant contribution directly stems from the dissipative nature of the system, and is present for all ${\bf k}$ in any parameter regime. We interpret this increase of the imaginary part as a dynamical instability which occurs even before the low-density solution becomes statically instable, e.g., before the density from Eq. (10) becomes imaginary. A stable gapless phase can be reached by increasing the coupling $J$ while keeping $j$ sufficiently weak, as illustrated in panel (c) of Fig. 2. In this case, the real part of the excitation frequencies vanishes on a circle around the M-point [${\bf k}a=(\pi,\pi)$], with linear dispersion relation around these gapless points. In summary, our analysis allows for distinguishing between four regions in parameter space: region of statical instability, region of dynamic instability, and a stable region divided into a gapped and a gapless regime (with respect to the real part of the excitation frequency).

From the Keldysh component of the Green function, we also obtain information about the momentum distribution and superfluid correlations. Not surprisingly, they directly reflect the different properties of the excitation spectra. As seen in Fig. 3, maxima in $n_{\bf k}$ and $|\Delta_{\bf k}|$ correspond to those points in ${\bf k}$-space where the real part of the excitation spectra exhibits minima, cf. Fig. 2. In the gapless phase, these maxima become sharply peaked. Divergencies in $n_{\bf k}$ and $|\Delta_{\bf k}|$ characterize the dynamically instable phase. Interestingly, upon tuning through a gapless point, the superfluid correlations $\Delta_{\bf k}$ exhibit a $\pi$-phase shift, see Fig. 3(f).

Transforming back to spatio-temporal variables, we further analyze the range of the equal-time correlations:

[TABLE]

The range of correlations is then quantified by the correlation length $\xi$ which we obtain by fitting the decay of correlations along a lattice direction to an exponential function:

[TABLE]

We have performed this analysis at $\gamma=0.2\mu$ and $U=0.5\mu$ in a weakly driven regime, $j=0.2\mu$, scanning several values of the cavity coupling $zJ/\mu$ in the interval $(0,2.5]$. According to our previous results, we expect quantitative and qualitative changes of the correlations near $zJ=\mu$. Indeed, the behavior can be described in the following way: For $zJ<\mu$, the spatial decay of is very well captured by the exponential fit, as seen in Fig. 4(a). When approaching $zJ\approx 1$ and for even larger values of $zJ$, the exponential decay is superposed by a spatial modulation of correlation, see Fig. 4(b–d). At the same time, the correlation length $\xi$ starts to increase, with ${\rm d}\xi/{\rm d}J$ being largest at $zJ/\mu=1$. This behavior is illustrated in Fig. 5(a). This observation supports the notion of two distinct regimes: At weak cavity coupling, $zJ<\mu$, the system exhibits only short-ranged correlations, while for $zJ>\mu$ the correlations acquire significantly longer range.

Finally, we have a look on the depletion. In order to justify the truncation of the action to second-order in $\delta$, we have to demand that the ratio between the density of flucuating particles, $n_{\delta}=\langle\delta_{0}\delta_{0}\rangle$, and the density of condensed particles $n=|b|^{2}$ is small. As seen in Fig. 5(b), at low drive $j=0.2\mu$, we have $n_{\delta}/n\lesssim 0.01$ for all values of $zJ/\mu$, rendering our description self-consistent. We also note that $n_{\delta}/n$ has a maximum near the critical coupling strength $zJ/\mu=1$ which might reflect a dynamical instability.

IV.2 High-density phase

In contrast to the low-density phase, the high-density solution to Eq. (10) does not appear to be further divided into different subphases. For any choice of parameters, we find gapped elementary excitations as illustrated in Fig. 6 (a). Again, the imaginary part of the excitation frequencies is always pinned to $i\gamma/2$. In contrast to the gapped phase at low density, with its lowest excitation at the M-point in the corner of the first Brillouin zone, the high density phase has its lowest excitation at the $\Gamma$-point in the center of the first Brillouin zone. This is also reflected by the position of the maximum of the momentum distribution $n_{\bf k}$ or superfluid correlations $\Delta_{\bf k}$, see Fig. 6 (b,c).

V Conclusions

We have studied a two-dimensional Bose-Hubbard model with dissipative losses and coherent driving. From the Schwinger-Keldysh action of the system, we obtain equations of motion which have one or several homogeneous solutions, leading to the distinction between mono- and multistable regimes. Fluctuations around these solutions yield the excitation spectra which allow for a further division of the low-density solution into a gapped, a gapless, and a dynamically instable regime. While our analysis was based on a mean-field approximation, the assumption of weak fluctuations was self-consistently confirmed. Our calculations may be relevant to future experiments with two-dimensional Bose-Hubbard models in a driven-dissipative environment. As an experimentally accessible figure of merit, we have computed the two-body correlations in the system. It is found that the transition from the gapped to the gapless regime is characterized through a significant increase of the correlation length.

Acknowledgements.

I acknowledge fruitful discussions with Bin Cao, Michael Foss-Feig, and Mohammad Maghrebi. Financial support from the NSF through PFC at JQI is acknowledged.