P.D. Nation and J.R. Johansson

For more information about QuTiP see http://qutip.org

In [1]:

```
import matplotlib.pyplot as plt
import numpy as np
from IPython.display import Image
from qutip import (about, destroy, hinton, ptrace, qdiags, qeye, steadystate,
tensor, wigner, wigner_cmap)
%matplotlib inline
```

In [2]:

```
Image(filename="images/optomechanical_setup.png", width=500, embed=True)
```

Out[2]:

\begin{equation}
\frac{\hat{H}}{\hbar}=-\Delta\hat{a}^{+}\hat{a}+\omega_{m}\hat{b}^{+}\hat{b}+g_{0}(\hat{b}+\hat{b}^{+})\hat{a}^{+}\hat{a}+E\left(\hat{a}+\hat{a}^{+}\right),
\end{equation}

The steady state density matrix of the optomechanical system plus the environment can be found from the Liouvillian superoperator $\mathcal{L}$ via

\begin{equation} \frac{d\rho}{dt}=\mathcal{L}\rho=0\rho, \end{equation}where $\mathcal{L}$ is typically given in Lindblad form \begin{align} \mathcal{L}[\hat{\rho}]=&-i[\hat{H},\hat{\rho}]+\kappa \mathcal{D}\left[\hat{a},\hat{\rho}\right]\\ &+\Gamma_{m}(1+n_{\rm{th}})\mathcal{D}[\hat{b},\hat{\rho}]+\Gamma_{m}n_{\rm th}\mathcal{D}[\hat{b}^{+},\hat{\rho}], \nonumber \end{align}

where $\Gamma_{m}$ is the coulping strength of the mechanical oscillator to its thermal environment with average occupation number $n_{th}$. As is customary, here we assume that the cavity mode is coupled to the vacuum.

Although, the steady state solution is nothing but an eigenvalue equation, the numerical solution to this equation is anything but trivial due to the non-Hermitian structure of $\mathcal{L}$ and its worsening condition number as the dimensionality of the truncated Hilbert space increases.

As of QuTiP version 3.0, the following steady state solvers are available:

**direct**: Direct LU factorization**eigen**: Calculates the eigenvector associated with the zero eigenvalue of $\mathcal{L}\rho$.**power**: Finds zero eigenvector using inverse-power method.**iterative-gmres**: Iterative solution via the GMRES solver.**iterative-lgmres**: Iterative solution via the LGMRES solver.**iterative-bicgstab**: Iterative solution via the BICGSTAB solver.**svd**: Solution via SVD decomposition (dense matrices only).

In [3]:

```
# System Parameters (in units of wm)
# -----------------------------------
Nc = 4 # Number of cavity states
Nm = 30 # Number of mech states
kappa = 0.3 # Cavity damping rate
E = 0.1 # Driving Amplitude
g0 = 2.4 * kappa # Coupling strength
Qm = 0.3 * 1e4 # Mech quality factor
gamma = 1 / Qm # Mech damping rate
n_th = 1 # Mech bath temperature
delta = -0.43 # Detuning
```

In [4]:

```
# Operators
# ----------
a = tensor(destroy(Nc), qeye(Nm))
b = tensor(qeye(Nc), destroy(Nm))
num_b = b.dag() * b
num_a = a.dag() * a
# Hamiltonian
# ------------
H = -delta * (num_a) + num_b + g0 * (b.dag() + b) * num_a + E * (a.dag() + a)
# Collapse operators
# -------------------
cc = np.sqrt(kappa) * a
cm = np.sqrt(gamma * (1.0 + n_th)) * b
cp = np.sqrt(gamma * n_th) * b.dag()
c_ops = [cc, cm, cp]
```

In [5]:

```
# all possible solvers
possible_solvers = ["direct", "eigen", "power", "iterative-gmres",
"iterative-bicgstab"]
# solvers used here
solvers = ["direct", "iterative-gmres"]
mech_dms = []
for ss in solvers:
if ss in ["iterative-gmres", "iterative-bicgstab"]:
use_rcm = True
else:
use_rcm = False
rho_ss, info = steadystate(
H,
c_ops,
method=ss,
use_precond=True,
use_rcm=use_rcm,
tol=1e-15,
return_info=True,
)
print(ss, "solution time =", info["solution_time"])
rho_mech = ptrace(rho_ss, 1)
mech_dms.append(rho_mech)
mech_dms = np.asarray(mech_dms)
```

direct solution time = 2.121961832046509 iterative-gmres solution time = 3.5488550662994385

In [6]:

```
for kk in range(len(mech_dms)):
c = np.where(
np.abs(mech_dms[kk].flatten() - mech_dms[0].flatten()) > 1e-5
)[0]
print("#NNZ for k = {}: {}".format(kk, len(c)))
```

#NNZ for k = 0: 0 #NNZ for k = 1: 0

`hinton()`

plot of the density matrix, we can see the magnitude of the diagonal elements is higher, such that the non-diagonal have a vanishing importance.

In [7]:

```
hinton(rho_mech, xlabels=[""] * Nm, ylabels=[""] * Nm);
```

`plt.spy()`

.

In [8]:

```
plt.spy(rho_mech.data, ms=1)
```

Out[8]:

<matplotlib.lines.Line2D at 0x7f9bc842e830>

In [9]:

```
diag = rho_mech.diag()
rho_mech2 = qdiags(diag, 0, dims=rho_mech.dims, shape=rho_mech.shape)
hinton(rho_mech2, xlabels=[""] * Nm, ylabels=[""] * Nm);
```

In [10]:

```
xvec = np.linspace(-20, 20, 256)
W = wigner(rho_mech2, xvec, xvec)
wmap = wigner_cmap(W, shift=-1e-5)
```

In [11]:

```
fig, ax = plt.subplots(figsize=(8, 6))
c = ax.contourf(xvec, xvec, W, 256, cmap=wmap)
ax.set_xlim([-10, 10])
ax.set_ylim([-10, 10])
plt.colorbar(c, ax=ax);
```

In [12]:

```
about()
```