Dissipative time crystal in an asymmetric nonlinear photonic dimer

We investigate the behavior of two coupled nonlinear photonic cavities, in the presence of inhomogeneous coherent driving and local dissipations. By solving numerically the quantum master equation, either by diagonalizing the Liouvillian superoperator or by using the approximated truncated Wigner approach, we extrapolate the properties of the system in a thermodynamic limit of large photon occupation. When the mean-field Gross-Pitaevskii equation predicts a unique parametrically unstable steady-state solution, the open quantum many-body system presents highly nonclassical properties and its dynamics exhibits the long-lived Josephson-like oscillations typical of dissipative time crystals, as indicated by the presence of purely imaginary eigenvalues in the spectrum of the Liouvillian superoperator in the thermodynamic limit.

Figure 1

(a) Phase diagram for the number and the nature of the GP steady-state solutions as a function of $J$ and $\tilde{F}$, for the fixed value of $\mathrm{\Delta }=2\kappa$. In the phase diagram, we distinguish the case of a single stable solution (1 sol. S), a single parametrically unstable solution (1 sol. P), three solutions (either with or without one parametrically unstable solution), and five solutions. The dashed line represents the value $J=1.2\kappa$, i.e., the value of $J$ considered in the results of Secs. 3 and 3. (b) and (c) Trajectories described by the mean fields ${\alpha }_1$ (red curve) and ${\alpha }_2$ (blue curve) according to the GP time evolution, for $\mathrm{\Delta }=2\kappa$, $J/\kappa =1.2$, and different values of $\tilde{F}$: panel (b) for $\tilde{F}=0.95$ (i.e., in the case of a single stable solution) and panel (c) for $\tilde{F}=1.5$ (i.e., in the case of a single parametrically unstable solution). The shown trajectories are obtained by numerically integrating the GP equations [Eq. (3)] up to $t\kappa ={10}^3$. The arrows indicate how the fields evolve for increasing time.

Figure 2

Region in the parameter space where GP predicts a single parametrically unstable solution. Each panel shows the results as a function of the parameters $J/\kappa$ and $\tilde{F}$, for different values of detuning: $\mathrm{\Delta }/\kappa$= 1 (a), 1.5 (b), 2 (c), 2.5 (d), 3 (e), and 3.5 (f).

Figure 3

Region in the parameter space where GP predicts a single parametrically unstable solution, for the case of a dimer made of resonators with different detuning frequencies [panels (a) and (b)], nonlinearities [panels (c) and (d)], or loss rates [panels (e) and (f)]. The coefficients for the first cavity are ${\mathrm{\Delta }}_1=2{\kappa }_1$, $U_1={\kappa }_1$, and $\tilde{F}=F\sqrt{U_1}/{\kappa }_1^{3/2}$. The coefficients for the second cavity are specified in each panel.

Figure 4

Rescaled steady-state expectation value of the photon occupations $n_1$ (a) and $n_2$ (b) in the two cavities versus the rescaled pump amplitude $\tilde{F}$. The lines with markers correspond to the different values of $U$. The solid, dashed, and dotted lines correspond to the predictions obtained from the steady-state solution of the GP equations, where the different line styles represent the different natures of the solution: stable (solid black line), one mode unstable (blue dotted line), and parametrically unstable (red dashed line). The light blue dot-dashed line represents the prediction for the photon occupation calculated from the time-averaged density matrix ${\widehat{\rho }}_{\mathrm{av}}$ [Eq. (6)]. The shaded area indicates approximatively the DTC phase. The other Hamiltonian parameters are $\mathrm{\Delta }/\kappa =2$ and $J/\kappa =1.2$.

Figure 5

Logarithmic negativity $E_N$ (a) and Von Neumann entropy $S$ (b) versus the rescaled pump amplitude and for different values of the nonlinearity $U$. The shaded area indicates approximatively the DTC phase. In panel (b), the light blue line represents the prediction for the von Neumann entropy calculated from the time-averaged density matrix ${\widehat{\rho }}_{\mathrm{av}}$ [Eq. (6)]. The other Hamiltonian parameters are $\mathrm{\Delta }/\kappa =2$ and $J/\kappa =1.2$.

Figure 6

Eigenvalues $\lambda$ of the Liouvillian superoperator, plotted in units of $\kappa$. Each panel shows the eigenvalues with the largest real part for a given pump amplitude $\tilde{F}$ and different values of nonlinearity $U$, corresponding to different colors and marker types. The black cross indicates the unique steady state. The different pump amplitudes corresponding to each panel are the following: $\tilde{F}=0.8$ (a), $\tilde{F}=1.0$ (b), $\tilde{F}=1.5$ (c), $\tilde{F}=1.8$ (d), $\tilde{F}=2.0$ (e), and $\tilde{F}=2.5$ (f). The other parameters are $\mathrm{\Delta }/\kappa =2$ and $J/\kappa =1.2$.

Figure 7

(a) Time evolution of the rescaled population difference $z$ between the two cavities for $\tilde{F}=1.5$, $\mathrm{\Delta }/\kappa =2$, $J/\kappa =1.2$, and different values of $U$. The results are obtained by averaging over ${10}^5$ stochastic trajectories obtained with the TWA. (b) Liouvillian gap $\mathrm{\Lambda }$ as a function of the nonlinearity $U$ for $\tilde{F}=1.5$, $\mathrm{\Delta }/\kappa =2$, and $J/\kappa =1.2$. The different symbols refer to different methods to extract the value of $\mathrm{\Lambda }$ [exact diagonalization or fit of TWA results for $z\left(t\right)$]. The dashed line represents a power-law fit of the data.

Figure 8

Time evolution of the population difference between the two cavities obtained for five different TWA trajectories, starting from the vacuum at $t=0$. The different panels show the results for different nonlinearities: $U=0.1$ (a), $U=0.05$ (b), and $U=0.01$ (c). The other parameters are $\tilde{F}=1.5$, $\mathrm{\Delta }/\kappa =2$, and $J/\kappa =1.2$.

Figure 9

Distribution of the TWA fields ${\alpha }_1$ (in red) and ${\alpha }_2$ (in blue) in phase space at different times $t$ and nonlinearities $U$. The black markers give the GP solution for the first ($+$) and second ($\times$) cavity for the given time. The distributions are obtained from the realization of $1.2\times {10}^4$ TWA trajectories.