Authors: J.R. Johansson and P.D. Nation

Modifications: C. Staufenbiel (2022)

In this tutorial we demonstrate the *Monte Carlo Solver* functionality implemented in `qutip.mcsolve()`

. For more information on the *MC Solver* refer to the QuTiP documentation.

We aim to reproduce the experimental results from:

Gleyzes et al., "Quantum jumps of light recording the birth and death of a photon in a cavity", Nature

446,297 (2007).

In particular, we will simulate the creation and annihilation of photons inside an optical cavity due to the thermal environment when the initial cavity is a single-photon Fock state $ |1\rangle$, as presented in Fig. 3 from the article.

First we import the relevant functionalities:

In [1]:

```
import matplotlib.pyplot as plt
import numpy as np
from qutip import about, basis, destroy, mcsolve, mesolve
%matplotlib inline
```

In this example, we consider a simple oscillator Hamiltonian $H = a^\dagger a$ and one initial photon in the cavity.

In [2]:

```
# number of modes in the cavity
N = 5
# Destroy operator
a = destroy(N)
# oscillator Hamiltonian
H = a.dag() * a
# Initial Fock state with one photon
psi0 = basis(N, 1)
```

The coupling to the external heat bath is described by a coupling constant $\kappa$ and the temperature of the heat bath is defined via the average photon number $\langle n \rangle$. In QuTiP the interaction between the system and heat bath is defined via the collapse operators. For this example, there are two collapse operators. One for photon annihilation ($C_1$) and one for photon creation ($C_2$):

$C_1 = \sqrt{\kappa (1 + \langle n \rangle)} \; a$

$C_2 = \sqrt{\kappa \langle n \rangle} \; a^\dagger$

We give some numerical values to the coupling constant $\kappa$ and the average photon number of the heat bath $\langle n \rangle$.

In [3]:

```
kappa = 1.0 / 0.129 # Coupling rate to heat bath
nth = 0.063 # Temperature with <n>=0.063
# collapse operators for the thermal bath
c_ops = []
c_ops.append(np.sqrt(kappa * (1 + nth)) * a)
c_ops.append(np.sqrt(kappa * nth) * a.dag())
```

The *Monte Carlo Solver* allows simulating an individual realization of the system dynamics. This is in contrast to e.g. the *Master Equation Solver*, which solves for the ensemble average over many identical realizations of the system. `qutip.mcsolve()`

also offers to average over many runs of identical system setups by passing the *number of trajectories* `ntraj`

to the function. If we choose `ntraj = 1`

the system is only simulated once and we see it's dynamics. If we choose a large value for `ntraj`

, the predictions will be averaged and therefore converge to the solution from `qutip.mesolve()`

.

We can also pass a list to `ntraj`

. `qutip.mcsolve()`

will calculate the results for the specified number of trajectories. Note that the entries need to be in ascending order, as the previous results are reused.

Here we are interested in the time evolution of $a^\dagger a$ for different numbers of `ntraj`

. We will compare the results to the predictions by `qutip.mesolve().

In [4]:

```
ntraj = [1, 5, 15, 904] # number of MC trajectories
tlist = np.linspace(0, 0.8, 100)
# Solve using MCSolve for different ntraj
mc = mcsolve(H, psi0, tlist, c_ops, [a.dag() * a], ntraj)
me = mesolve(H, psi0, tlist, c_ops, [a.dag() * a])
```

Using the above results we can reproduce Fig. 3 from the article mentioned above. The individual figures plot the time evolution of $\langle a^\dagger a \rangle$ for the system we set up above. The effect of using different `ntraj`

for the simulation using `mcsolve`

is shown. When choosing `ntraj = 1`

we see the dynamics of one particular quantum system. If `ntraj > 1`

the output is averaged over the number of realizations.

In [5]:

```
fig = plt.figure(figsize=(8, 8), frameon=False)
plt.subplots_adjust(hspace=0.0)
for i in range(4):
ax = plt.subplot(4, 1, i + 1)
ax.plot(
tlist, mc.expect[i][0], "b", lw=2,
label="#trajectories={}".format(ntraj[i])
)
ax.plot(tlist, me.expect[0], "r--", lw=2)
ax.set_yticks([0, 0.5, 1])
ax.set_ylim([-0.1, 1.1])
ax.set_ylabel(r"$\langle P_{1}(t)\rangle$")
ax.legend()
ax.set_xlabel(r"Time (s)");
```

In [6]:

```
about()
```

In [7]:

```
np.testing.assert_allclose(me.expect[0], mc.expect[3][0], atol=10**-1)
assert np.all(np.diff(me.expect[0]) <= 0)
```