Iterative solutions to the steady-state density matrix for optomechanical systems

We present a sparse matrix permutation from graph theory that gives stable incomplete lower-upper preconditioners necessary for iterative solutions to the steady-state density matrix for quantum optomechanical systems. This reordering is efficient, adding little overhead to the computation, and results in a marked reduction in both memory and runtime requirements compared to other solution methods, with performance gains increasing with system size. Either of these benchmarks can be tuned via the preconditioner accuracy and solution tolerance. This reordering optimizes the condition number of the approximate inverse and is the only method found to be stable at large Hilbert space dimensions. This allows for steady-state solutions to otherwise intractable quantum optomechanical systems.

Figure 1

(a) Sparse matrix structure for the modified Liouvillian $\tilde{\mathcal{L}}$ of an optomechanical system with $N_c=4$ and $N_m=8$ along with the matrix bandwidth $B$ and profile $P$. Here it is assumed that the mechanical oscillator is in a nonzero thermal environment. The elements corresponding to matrix $\mathcal{T}$ (red), here enlarged for visibility, responsible for the large total bandwidth, are along the upper row. (b) RCM reordering ${\tilde{\mathcal{L}}}_{\mathrm{RCM}}$, where the permuted matrix has a total bandwidth $\sim 5$ times lower than the naturally ordered matrix, and a profile reduction of $\sim 2.5$. (c) Nonsymmetric COLAMD ordering of $\tilde{\mathcal{L}}$. Only the structure of $\tilde{\mathcal{L}}$ is required for these reorderings. The number of nonzero elements is $N_{\mathrm{NZ}}\left(\tilde{\mathcal{L}}\right)=8957$.

Figure 2

(a) Log-plot of the bandwidth reduction factor after RCM reordering for an optomechanical system with mechanical oscillator coupled to a nonzero thermal environment for a varying number of cavity and mechanical resonator states, $N_c$ and $N_m$, respectively. (b) Log-plot of the corresponding profile reduction factor.

Figure 3

Fill-factors for the direct LU decomposition of the modified Liouvillian $\tilde{\mathcal{L}}$ using natural (gray), RCM (red), and COLAMD (blue) matrix ordering as a function of Hilbert space dimensionality. Here the system parameters (in units of ${\omega }_m$) are $\kappa =0.05$, $\Delta =-\kappa$, $g_0=3\kappa$, $E=0.25$, $Q_m={10}^4$, and $n_{\mathrm{th}}=31$.

Figure 4

(a) Fill-factors as a function of detuning for direct LU factorization using RCM (brown) or COLAMD (gray) reordering, along with fill-factors for iLU preconditioners based on ${\tilde{\mathcal{L}}}_{\mathrm{RCM}}$ at $n_{\mathrm{th}}=0$ (blue), $n_{\mathrm{th}}={10}^{-15}$ (green), and $n_{\mathrm{th}}=3$ (red). Those based on the COLAMD ordering of $\tilde{\mathcal{L}}$ at $n_{\mathrm{th}}=0$ (yellow) and $n_{\mathrm{th}}=3$ (purple) are also presented. (b) Number of iterations needed to reach machine precision for the iLU preconditioners constructed in (a). Only those points where the solution converged within 1000 iterations are displayed. (c) Log-plot of the estimated condition number $\left|\right|\mathcal{M}\vec{e}{\left|\right|}_{\infty }$ for the iLU preconditioners given in (a). (d) Computational speedup of iterative solvers as measured against the runtime of direct LU factorization using COLAMD reordering. Only those points where convergence was reached are displayed. For (a)–(d), the number of cavity and mechanical oscillator states is $N_c=4$ and $N_m=200$, while the system parameters (in units of ${\omega }_m$) are $\kappa =0.2$, $g_0=2.5\kappa$, $Q_m={10}^4$, and $E=0.1$. The iLU drop tolerance and allowed fill-in are $d={10}^{-4}$ and $p=100$, respectively. [(e)–(h)] Results of iterative simulations for the system parameters from Ref. [13], and given in Fig. 3, for $N_c=10$ and $N_m=50$. Here the iLUTP parameters are set to $d=5\times {10}^{-5}$ and $p=300$. Only the direct LU (brown) and iterative (red) RCM results for $n_{\mathrm{th}}=31$, along with the corresponding COLAMD direct (gray) and iterative (purple) simulations, are presented.

Figure 5

(a) Wigner distribution for the mechanical oscillator density matrix at (in units of ${\omega }_m$): $E=0.45$, $\kappa =0.05$, $\Delta =-2\kappa$, $g_0=0.32$, and $n_{\mathrm{th}}=20$. Here $x_{\mathrm{zp}}$ and $p_{\mathrm{zp}}$ are the zero-point amplitudes of the oscillator's position and momentum, respectively. (b) Phonon number distribution (gray) for the oscillator given in (a) along with fitted coherent states for the lower (red) and upper (blue) limit cycles. The total Fano factor for the distribution is also given. Simulation parameters are $N_c=5$, $N_m=160$, $d={10}^{-4}$, and $p=200$. Only the first 125 oscillator Fock states are displayed for better limit-cycle visibility.