Computation of the asymptotic states of modulated open quantum systems with a numerically exact realization of the quantum trajectory method

Quantum systems out of equilibrium are presently a subject of active research, both in theoretical and experimental domains. In this work, we consider time-periodically modulated quantum systems that are in contact with a stationary environment. Within the framework of a quantum master equation, the asymptotic states of such systems are described by time-periodic density operators. Resolution of these operators constitutes a nontrivial computational task. Approaches based on spectral and iterative methods are restricted to systems with the dimension of the hosting Hilbert space $\dim \mathcal{H}=N\lesssim 300$, while the direct long-time numerical integration of the master equation becomes increasingly problematic for $N\gtrsim 400$, especially when the coupling to the environment is weak. To go beyond this limit, we use the quantum trajectory method, which unravels the master equation for the density operator into a set of stochastic processes for wave functions. The asymptotic density matrix is calculated by performing a statistical sampling over the ensemble of quantum trajectories, preceded by a long transient propagation. We follow the ideology of event-driven programming and construct a new algorithmic realization of the method. The algorithm is computationally efficient, allowing for long “leaps” forward in time. It is also numerically exact, in the sense that, being given the list of uniformly distributed (on the unit interval) random numbers, $\{{\eta }_1,{\eta }_2,...,{\eta }_n\}$, one could propagate a quantum trajectory (with ${\eta }_i$'s as norm thresholds) in a numerically exact way. By using a scalable $N$-particle quantum model, we demonstrate that the algorithm allows us to resolve the asymptotic density operator of the model system with $N=2000$ states on a regular-size computer cluster, thus reaching the scale on which numerical studies of modulated Hamiltonian systems are currently performed.

Figure 1

Floatlike performance of Algorithm 2 during the propagation of a quantum trajectory. Every local maximum in the dependence $\delta t$ vs. number of steps (minimum in the depth $s$ dependence) indicates an occurrence of a jump after which the algorithm performs a chain of bisections to reach the maximal depth $S=20$. After every step during which no jump occurred, the algorithm doubles the step size. The average time between two consecutive jumps (red line) is the average height of the local maximum minus 1. The two sequences were monitored during the sampling of the asymptotic state.

Figure 2

(a) Structure of the sampled stroboscopic density matrix ${\varrho }^{\mathrm{att}}\left(0\right)$ ($N=1024$) and (b) time evolution of the mean ${\overline{\varrho }}_{11}\left(t\right)$ (thick blue line) and variance $var\left[{\varrho }_{11}\left(t\right)\right]$ (thin red line) for $t\in [40T,50T]$. The system size is $N=256$ and the sampling was performed over ${10}^5$ independent trajectories initiated at the state $|{\psi }^{\mathrm{init}}\rangle =|1\rangle$ and then propagated to the time $t_p=50T$. The inset depicts the limit-cycle evolution of the means and variances for two diagonal elements, ${\varrho }_{1,1}$ and ${\varrho }_{128,128}$, during one period of modulations, $t\in [1000T,1001T]$. Curves for later periods are indistinguishable from the presented ones. The parameters are the same as in Fig. 1.

Figure 3

Spectral norm of the difference between the density matrix stroboscopically sampled with quantum trajectory algorithm and numerically exact asymptotic solution, $\epsilon =‖{\overline{\varrho }}^{\mathrm{att}}\left(mT\right)-{\varrho }^{ex}(\tau =0)‖$, for the two-mode model, Eqs. (6) and (7). Here $N=100$, $J=1$, ${\mu }_0=1.5$, ${\mu }_1=1$, $U=3$, ${\gamma }_0=0.1$.

Figure 4

Attractors of the two-mode model, Eqs. (6) and (7), with $N-1=1023$ bosons. Left panel: Husimi distribution $p(\vartheta ,\phi )$ of the stroboscopic density matrix ${\varrho }^{\mathrm{att}}\left(0\right)$ (top) and Poincaré plot of the attractor of the corresponding mean-field system, Eq. (8) (bottom). The density matrix was sampled with ${10}^5$ stroboscopic realizations. Right panel: Plane version of $p(\vartheta ,\phi )$ (top) compared to the probability density function (pdf) $P(\vartheta ,\phi )$ of the classical mean-field attractor (bottom). The pdf was estimated by sampling the histogram on a $200\times 200$ bin grid with ${10}^8$ stroboscopic points. Note that, because of a high localization of the density near point $(0,3\pi /2)$, the logarithmic scale is used. The parameters are $J=1$, ${\mu }_0=1.5$, ${\mu }_1=1$, $U=3$, ${\gamma }_0=0.1$.