Efficient choice of colored noise in the stochastic dynamics of open quantum systems

The stochastic Liouville–von Neumann (SLN) equation describes the dynamics of an open quantum system reduced density matrix coupled to a non-Markovian harmonic environment. The interaction with the environment is represented by complex colored noises which drive the system, and whose correlation functions are set by the properties of the environment. We present a number of schemes capable of generating colored noises of this kind that are built on a noise amplitude reduction procedure [Imai et al., Chem. Phys. 446, 134 (2015)], including two analytically optimized schemes. In doing so, we pay close attention to the properties of the correlation functions in Fourier space, which we derive in full. For some schemes the method of Wiener filtering for deconvolutions leads to the realization that weakening causality in one of the noise correlation functions improves numerical convergence considerably, allowing us to introduce a well-controlled method for doing so. We compare the ability of these schemes, along with an alternative optimized scheme [Schmitz and Stockburger, Eur. Phys. J.: Spec. Top. 227, 1929 (2019)], to reduce the growth in the mean and variance of the trace of the reduced density matrix, and their ability to extend the region in which the dynamics is stable and well converged for a range of temperatures. By numerically optimizing an additional noise scaling freedom, we identify the scheme which performs best for the parameters used, improving convergence by orders of magnitude and increasing the time accessible by simulation.

Figure 1

Average magnitude of $\nu \left(t\right)$ taken across 500 realizations for each $t_{max}$ for (a) the constrained and (b) reduced schemes with (green) and without (blue) the Wiener filter using $\gamma =0.01$, $\beta =1$, $\mathrm{\Delta }=1$, $\epsilon =-1$, $\alpha =0.05$, and ${\omega }_c=25$.

Figure 2

The $\eta \text{−}\nu$ correlation function with different values of the parameter $\gamma$ in the Wiener filter for the constrained noise scheme. The $\eta \nu$ optimized scheme overlaps with the desired correlation $K_{\eta \nu }$ such that $K_{\eta \nu }$ could not be seen, so it is not shown. $\beta =1$, $t_{max}=12$, $dt=0.01$, ${\omega }_c=25$, $\alpha =0.05$ for ${10}^4$ realizations. The zoomed inset highlights the region in which the symmetrization of the correlation as $\gamma$ increases can be clearly seen.

Figure 3

Comparison of the calculated expectation of the $z$-spin $\langle {\sigma }_z\left(t\right)\rangle$ using the constrained scheme with a range of $\gamma$ values to control deconvolution, and the $\eta \nu -$optimized scheme $\langle {\sigma }_z^{\eta \nu }\left(t\right)\rangle$, which minizes the combined magnitude of $\eta$ and $\nu$. (a) The system is initialized with the $z$-spin being 1 and all other spins being zero and is evolved in time with a constant Hamiltonian, so that the $z$-spin relaxes to its equilibrium value. The canonical equilibrium value of 0.05 for the spin, obtained using the imaginary time evolution methods outlined in Ref. [18], is shown (dashed line) to confirm the validity of the optimized scheme. The inset shows the deviation of the $z$-spin for the constrained scheme from the optimized one. A zoomed-in area at the final stages of equilibration shows the $z$-spin in detail. (b) A Landau-Zener sweep with a time-dependent Hamiltonian where the system is driven linearly with $\epsilon \left(t\right)=5t$, approaching a known asymptotic limit (dashed line). The system is initialized at $t=-5$ with the $z$-spin equal to 1 and all other spins being zero, using the modified Landau-Zener limit of 0.516 to account for the finiteness of the simulation window as outlined in [18]. For clarity, data are no longer plotted once they exceed the vertical scales shown, as the solution becomes unstable and grows exponentially. $\beta =0.1$, $dt={10}^{-2}$, $\mathrm{\Delta }=1$, $\epsilon =-1$, $\alpha =0.05$, and ${\omega }_c=25$ for ${10}^6$ realizations.

Figure 4

The total absolute deviation of the dynamics produced using the constrained scheme for different $\gamma$ from the dynamics produced using the $\eta \nu$-optimized scheme. Values were calculated using the data shown in the insets of Figs. 3 and 3. Results for $\gamma =0$ and 0.001 are not shown because the dynamics is diverging so the deviation is very large. The value of $\gamma$ which minimizes the total deviation for both the constant Hamiltonian [open circles, Fig. 3] and the Landau-Zener sweep [filled circles, Fig. 3] is 0.01.

Figure 5

Dynamics of the (01) element of the reduced density matrix according to the quantum nondemolition model, showing the real part $Re\left[\langle {\rho }_{01}\left(t\right)\rangle \right]$ in (a), and the imaginary part $\mathrm{Im}\left[\langle {\rho }_{01}\left(t\right)\rangle \right]$ in (b). The exact solution (black line) is compared to the SLN numerical solutions for 1000 (blue line) and 50 000 (green line) realizations, using the $\eta \nu$-optimized scheme with optimal scaling $\lambda =0.5$ (see Sec. 4c).

Figure 6

(a)–(c) The mean value of the trace, $\langle Tr(\mathit{\rho }\left(t)\right)\rangle$, calculated using different schemes with the system having been initialized in the state $|1\rangle$ for a constant Hamiltonian. The absolute value of $\langle Tr(\mathit{\rho }\left(t\right))\rangle$ is shown so that the linear growth on a logarithmic scale is clear. The insets highlight the timescales on which the simulation is numerically stable, showing $\langle Tr(\mathit{\rho }\left(t\right))\rangle$ directly; schemes are not shown for timescales beyond which the trace is clearly diverging. (d–f) The variance of the trace and (g)–(i) the standard error of the mean trace calculated over time windows which were 100 time steps long. For columns read from left to right, the inverse temperature increases as $\beta =0.1,1,10$, respectively. The like (red), constrained (blue), and reduced (green) schemes are shown, as well as the $\nu$-optimized scheme which minimizes $\langle |\nu \left(t\right){|}^2\rangle$ (black solid), $\eta \nu$-optimized scheme which minimizes the sum ${\langle \left|\eta \left(t\right)\right|}^2{\rangle +\langle \left|\nu \left(t\right)\right|}^2\rangle$ (black dashed), and the application of the convex optimization scheme [21] as implemented using Eqs. (37, 38, 39, 40, 41, 42) (yellow). All calculations have been done using the same system as in Fig. 3. $\mathrm{\Delta }=1$, $\epsilon =-1$, $\alpha =0.05$, $\mathrm{\Delta }t={10}^{-2}$, ${\omega }_c=25$, and ${10}^5$ realizations. No rescaling of the noises was employed.

Figure 7

(a) The standard error of $\langle Tr\mathit{\rho }\left(t\right)\rangle$ at its final time step $t_{max}$ as a function of the scaling factor $\lambda$ for several values of inverse temperature $\beta$ and coupling strength $\alpha$. For each scaling factor, 1000 runs for real time dynamics were performed. $t_{max}$ = 40, $dt={10}^{-3}$, ${\omega }_c=25$, and $\mathrm{\Delta }=1$, $\epsilon =-1$ for $\alpha =0.05,0.1$ and $\beta =0.1,1,10.$ In this case, the optimum value of $\lambda$ which minimizes the growth of the trace is $\approx 0.5$. (b) The variance of the trace having used the $\eta \nu$-optimized scheme with scaling, with a desired ratio between $f_2\left(t\right)$ and $g_2\left(t\right)$ of $\lambda =0.5$ (solid lines) for $\beta =0.1,1,10$, with the mean trace shown in the inset. The convex optimized scheme (dashed lines) has been reproduced here from Figs. 6 for comparison.