Stabilizing strongly correlated photon fluids with non-Markovian reservoirs

We introduce a frequency-dependent incoherent pump scheme with a square-shaped spectrum as a way to study strongly correlated photons in arrays of coupled nonlinear resonators. This scheme can be implemented via a reservoir of population-inverted two-level emitters with a broad distribution of transition frequencies. Our proposal is predicted to stabilize a nonequilibrium steady state sharing important features with a zero-temperature equilibrium state with a tunable chemical potential. We confirm the efficiency of our proposal for the Bose-Hubbard model by computing numerically the steady state for finite system sizes: first, we predict the occurrence of a sequence of incompressible Mott-insulator-like states with arbitrary integer densities presenting strong robustness against tunneling and losses. Secondly, for stronger tunneling amplitudes or noninteger densities, the system enters a coherent regime analogous to the superfluid state. In addition to an overall agreement with the zero-temperature equilibrium state, exotic nonequilibrium processes leading to a finite entropy generation are pointed out in specific regions of parameter space. The equilibrium ground state is shown to be recovered by adding frequency-dependent losses. The promise of this improved scheme in view of quantum simulation of the zero-temperature many-body physics is highlighted.

Figure 1

Panel (a): plot for various values of ${\mathrm{\Delta }}_{\mathrm{em}}$ of the “square-shaped” emission spectrum ${\mathcal{S}}_{\mathrm{em}}\left(\omega \right)/{\mathrm{\Gamma }}_{\mathrm{em}}^0=s_{\mathrm{em}}\left(\omega \right)$ [defined in Eq. (9)] in units of ${\mathrm{\Gamma }}_{\mathrm{em}}^0$. Panel (b): average photon number $n_{\mathrm{ph}}$ as a function of $\mu ={\omega }_{+}-{\omega }_{\mathrm{cav}}$ (i.e., varying ${\omega }_{\mathrm{cav}}$) for a single-site system with various values of ${\mathrm{\Delta }}_{\mathrm{em}}$. Parameters of panel (b): ${\mathrm{\Gamma }}_{\mathrm{em}}^0/U=3\times {10}^{-4}, {\mathrm{\Gamma }}_{\mathrm{l}}/U={10}^{-5}$, and ${\omega }_{-}/U=-40$.

Figure 2

Steady-state properties for a limited choice of parameters. Panels (a),(c),(e) [panels (b),(d),(f)]: average steady-state photon number per site $n_{\mathrm{ph}}$ (condensed fraction $x_{\mathrm{BEC}}$), for $L=2,5,7$, respectively. Panel (g) [panel (h)]: average $n_{\mathrm{ph}}$ ($x_{\mathrm{BEC}}$) in a $T=0$ equilibrium system, for $L=7$. Panel (i) [panel (j)]: steady-state particle number relative fluctuations $\mathrm{\Delta }n$ [entropy $S=-\langle \text{ln}\left({\rho }_{\infty }\right)\rangle$], for $L=7$. In panels (g) and (j) dash-dotted black lines indicate the MPS $T=0$ prediction for the first Mott lobes. Parameters used in all panels [except (g),(h)]: ${\mathrm{\Gamma }}_{\mathrm{l}}/{\mathrm{\Gamma }}_{\mathrm{em}}^0={10}^{-3}, U/{\mathrm{\Delta }}_{\mathrm{em}}={10}^6, {\mathrm{\Gamma }}_{\mathrm{em}}^0/{\mathrm{\Delta }}_{\mathrm{em}}={10}^{-2}, {\omega }_{+}=0$, and ${\omega }_{-}/{\mathrm{\Delta }}_{\mathrm{em}}=-4\times {10}^7$. Cutoff in particle number per site $N_{\max }=6$ [panels (a),(b)] and $N_{\max }=3$ [panels (c)–(j)].

Figure 3

Steady-state properties for $L=3$, for state-of-the-art parameters in circuit QED. Panel (a): average photon number per site $n_{\mathrm{ph}}$. Panel (b): condensed fraction $x_{\mathrm{BEC}}$. Panel (c): relative fluctuations of the total particle number $\mathrm{\Delta }n$. Parameters inspired from circuit-QED systems [35, 44]: $U=200\times 2\pi \text{MHz}, {\mathrm{\Delta }}_{\mathrm{em}}=0.5\times 2\pi \text{MHz}, {\mathrm{\Gamma }}_{\mathrm{em}}^0=30\times 2\pi \text{kHz}$, and ${\mathrm{\Gamma }}_{\mathrm{l}}=1\times 2\pi \text{kHz}$. In order to be able to correctly see the higher lobes, we had to increase the maximum allowed number of particles per site to $N_{\max }=4$, and correspondingly to reduce the system size to $L=3$.

Figure 4

Steady-state statistical properties for $L=5$ at fixed $\mu =0.55U$. Panel (a): fidelity $\mathcal{F}$ between the steady state and the Hamiltonian ground state (blue solid line), occupancy ${\pi }_0$ of the most populated state ${|\psi \rangle }_{+}$ (purple crosses), overlap $|\langle {\psi }_{+}|\mathrm{GS}\rangle {|}^2$ between the ground state and the most populated state (red dotted line), and entropy $S$ (green dashed line). Panel (b): number fluctuations $\mathrm{\Delta }n$ (green solid line) and condensed fraction (blue circles), compared to the $T=0$ equilibrium value (orange dash-dotted line). Same parameters as in Fig. 2. In order to be able to perform exact diagonalization of the Liouvillian and avoid using the secular approximation, we had to choose a smaller system size $L=5$ as compared to Fig. 2.

Figure 5

Left panel: occupancy ${\pi }_0$ of the most populated quantum state at steady state for a five sites system. Central (right) panel: spectrum of $H_{\text{eff}}$ at the values $J/U=0.1$ and $\mu /U=0.38$ ($\mu /U=0.55$) indicated by the point $\left[\alpha \right]$ ($\left[\beta \right]$) in the left panel. Parameters used: ${\mathrm{\Gamma }}_{\mathrm{l}}/{\mathrm{\Gamma }}_{\mathrm{em}}={10}^{-3}, U/{\mathrm{\Delta }}_{\text{em}}={10}^6, {\mathrm{\Gamma }}_{\mathrm{em}}^0/{\mathrm{\Delta }}_{\text{em}}={10}^{-2}, {\omega }_{+}=0$, and ${\omega }_{-}/{\mathrm{\Delta }}_{\text{em}}=-4\times {10}^7$.

Figure 6

Steady-state properties using frequency-dependent losses with a square spectrum in addition to frequency-dependent emission. Simulations were done for a $L=7$ sites periodic chain. Panel (a) [panel (c)]: average steady-state photon number per site $n_{\mathrm{ph}}$ (condensed fraction $x_{\mathrm{BEC}}$). Panel (b): steady-state entropy. Panel (d): $1-\mathcal{F}$, where $\mathcal{F}=\langle \mathrm{GS}\left|{\rho }_{\infty }\right|\mathrm{GS}\rangle$ is the fidelity of the steady-state density matrix ${\rho }_{\infty }$ with the ground state $|\mathrm{GS}\rangle$ of the Hamiltonian $H_{\mathrm{eff}}$, i.e., a $T=0$ state of chemical potential $\mu$. Same parameters as in Fig. 2 except for the additional frequency-dependent losses: ${\mathrm{\Gamma }}_{\mathrm{L}}^0={\mathrm{\Gamma }}_{\mathrm{em}}^0, {\omega }_{\mathrm{L}}-{\omega }_{\mathrm{cav}}={\omega }_{+}-{\omega }_{-}$, and ${\mathrm{\Delta }}_{\mathrm{L}}={\mathrm{\Delta }}_{\mathrm{em}}$.

Figure 7

Steady-state properties using only one or two emitting sites in different geometries. Panel (a) [(c)] shows the steady-state average density $n_{\mathrm{ph}}$, and panel (b) [(d)] the entropy for a $L=4$ sites system with periodic boundary conditions with emitters localized on the first two sites (an open chain with emitters only on the first site). Same parameters as in Fig. 2 except for the stronger emission rate to compensate the reduced number of emitting sites: ${\mathrm{\Gamma }}_{\mathrm{em}}^0/{\mathrm{\Delta }}_{\mathrm{em}}=\frac{L}{2}\times {10}^{-2}$ in panels (a),(b) [$L\times {10}^{-2}$ in panels (c),(d)]. A smaller lattice of $L=4$ sites had to be used because of the broken translational invariance.

Figure 8

Temporal profile of the time-dependent emitter transition frequency required to mimic the square-shaped emission spectrum.

Figure 9

Time evolution of a single-cavity configuration using a single temporally modulated emitter to mimic the effect of a square emission spectrum. Different panels from left to right refer to increasing values of the modulation speed $v_{\omega }=\left|\frac{d{\omega }_{at}}{dt}\right|$. Upper panels: average photon number $〈N〉\left(t\right)$. Lower panels: relative fluctuations of the total particle number $\mathrm{\Delta }n\left(t\right)=\sqrt{〈N^2〉\left(t\right)-{〈N〉}^2\left(t\right)}/〈N〉\left(t\right)$. The blue lines ending at long times at $n\simeq 1$ and higher $\mathrm{\Delta }n$ (red lines ending at $n\simeq 2$ and lower $\mathrm{\Delta }n$) correspond to a single-photon (two-photon) Fock state stabilization. The photon-emitter Rabi coupling ${\mathrm{\Omega }}_R$ was chosen for each modulation speed in such a way as to set the stationary value of the photon number close to the desired occupation number: $〈N〉\left(t\right){\simeq }_{t\to \infty }1$ (2). Parameters inspired from state-of-the-art in circuit-QED systems [35, 44]: $U=200\times 2\pi \text{MHz}, {\mathrm{\Gamma }}_{\mathrm{p}}=0.5\times 2\pi \text{MHz}$, and ${\mathrm{\Gamma }}_{\mathrm{l}}=1\times 2\pi \text{kHz}$. For the single-photon (two-photon) simulation $\mu ={\omega }_{+}-{\omega }_{\mathrm{cav}}=U/2$ ($3U/2)$ and ${\omega }_{+}-{\omega }_{-}=0.6U$ ($1.6U$). Choice for the modulation speed, from left to right: $v_{\omega }=7.5\times {\left(2\pi \text{MHz}\right)}^2, 15\times {\left(2\pi \text{MHz}\right)}^2, 30\times {\left(2\pi \text{MHz}\right)}^2$, and $50\times {\left(2\pi \text{MHz}\right)}^2$. Correspondingly, for the single-photon simulation: ${\mathrm{\Omega }}_R=0.83\times 2\pi \text{Hz}, 0.9\times 2\pi \text{Hz}, 0.97\times 2\pi \text{Hz}$, and $1.07\times 2\pi \text{Hz}$ (for the two-photon simulation: ${\mathrm{\Omega }}_R=1\times 2\pi \text{Hz}, 1.07\times 2\pi \text{Hz}, 1.15\times 2\pi \text{Hz}$, and $1.24\times 2\pi \text{Hz}$).