Andrianov et al. [1] resort to the quantum Monte Carlo (QMC) simulations to their Lindblad Eq. (1) due to the “no man’s land” regime of model parameters. However, the Liouvillian superoperator of their Eq. (1) has sparse representation in the Fock space, regardless of any parameter values, and thus can be targeted by standard sparse numerical methods implemented in most numerical packages, including Matlab, Mathematica, or QuTiP [2]. We employed direct calculations of the stationary state as well as the photonic mode spectrum for the parameters ${\omega }_a={\omega }_b=0,{\mathrm{\Omega }}_R=0.35,{\gamma }_a=0.002,{\gamma }_b=0.25,{\mathrm{\Omega }}_b=9{\mathrm{\Omega }}_a$ and a range of ${\mathrm{\Omega }}_a$’s used in the original article (Private communication with N. M. Chtchelkatchev; these parameters cannot be figured out from the original text due to omissions and misprints.). These calculations take several seconds on a standard laptop for more-than-sufficient photon numbers up to 250. We recovered the original results for N and Q in Fig. 2 with the identification $\mathrm{\Omega }\equiv 2{\mathrm{\Omega }}_a$.

However, we cannot confirm the plots for the Wigner functions of Figs. 3(g)–3(j) and the non-monotonic behavior of the spectral line for the large drives in Fig. 4. The corresponding results are shown in our Figs. 1(a)–1(d) (Wigner function) and Figs. 1(e)–1(h) and 1(i) (spectra). They do not exhibit any fading of one of the Wigner function peaks or broadening of the spectral line with increasing drive. Also, the absolute value of the half width comes out wrong even for weaker drives. Our results are confirmed by several methods including (numerically inefficient) quantum trajectory simulations in QuTiP and strongly suggest serious convergence issues in the original QMC simulations. This is consistent with the calculated first nonzero eigenvalue of the Liouvillian ${\lambda }_1$ having a large spectral gap from the rest and determining the spectral line half width, cf., Fig. 1(i). The existence of an emergent long timescale $1/{\lambda }_1$, which is a strong hint for bistability [3], may explain the problems with the QMC convergence. These spotted errors challenge the proposed physical interpretation of the model’s behavior.

 

Fig. 1. Wigner function (a)–(d), spectrum (e)–(h), and line half width (i) corresponding to Figs. 3 and 4 in Ref. [1].

Download Full Size | PDF