# Monte Carlo for Non-Markovian Dynamics

The Monte Carlo solver of QuTiP can also be used to solve the dynamics of time-local non-Markovian master equations, i.e., master equations of the Lindblad form

(1)\[\dot\rho(t) = -\frac{i}{\hbar} [H, \rho(t)] + \sum_n \frac{\gamma_n(t)}{2} \left[2 A_n \rho(t) A_n^\dagger - \rho(t) A_n^\dagger A_n - A_n^\dagger A_n \rho(t)\right]\]

with “rates” \(\gamma_n(t)\) that can take negative values. This can be done with the `nm_mcsolve` function. The function is based on the influence martingale formalism [Donvil22] and formally requires that the collapse operators \(A_n\) satisfy a completeness relation of the form

(2)\[\sum_n A_n^\dagger A_n = \alpha \mathbb{I} ,\]

where \(\mathbb{I}\) is the identity operator on the system Hilbert space and \(\alpha>0\). Note that when the collapse operators of a model don’t satisfy such a relation, `nm_mcsolve` automatically adds an extra collapse operator such that (2) is satisfied. The rate corresponding to this extra collapse operator is set to zero.

Technically, the influence martingale formalism works as follows. We introduce an influence martingale \(\mu(t)\), which follows the evolution of the system state. When no jump happens, it evolves as

(3)\[\mu(t) = \exp\left( \alpha\int_0^t K(\tau) d\tau \right)\]

where \(K(t)\) is for now an arbitrary function. When a jump corresponding to the collapse operator \(A_n\) happens, the influence martingale becomes

(4)\[\mu(t+\delta t) = \mu(t)\left(\frac{K(t)-\gamma_n(t)}{\gamma_n(t)}\right)\]

Assuming that the state \(\bar\rho(t)\) computed by the Monte Carlo average

(5)\[\bar\rho(t) = \frac{1}{N}\sum_{l=1}^N |\psi_l(t)\rangle\langle \psi_l(t)|\]

solves a Lindblad master equation with collapse operators \(A_n\) and rates \(\Gamma_n(t)\), the state \(\rho(t)\) defined by

(6)\[\rho(t) = \frac{1}{N}\sum_{l=1}^N \mu_l(t) |\psi_l(t)\rangle\langle \psi_l(t)|\]

solves a Lindblad master equation with collapse operators \(A_n\) and shifted rates \(\gamma_n(t)-K(t)\). Thus, while \(\Gamma_n(t) \geq 0\), the new “rates” \(\gamma_n(t) = \Gamma_n(t) - K(t)\) satisfy no positivity requirement.

The input of `nm_mcsolve` is almost the same as for `mcsolve`. The only difference is how the collapse operators and rate functions should be defined. `nm_mcsolve` requires collapse operators \(A_n\) and target “rates” \(\gamma_n\) (which are allowed to take negative values) to be given in list form `[[C_1, gamma_1], [C_2, gamma_2], ...]`. Note that we give the actual rate and not its square root, and that `nm_mcsolve` automatically computes associated jump rates \(\Gamma_n(t)\geq0\) appropriate for simulation.

We conclude with a simple example demonstrating the usage of the `nm_mcsolve` function. For more elaborate, physically motivated examples, we refer to the accompanying tutorial notebook. Note that the example also demonstrates the usage of the `improved_sampling` option (which is explained in the guide for the Monte Carlo Solver) in `nm_mcsolve`.

    times = np.linspace(0, 1, 201)
    psi0 = basis(2, 1)
    a0 = destroy(2)
    H = a0.dag() * a0

    # Rate functions
    gamma1 = "kappa * nth"
    gamma2 = "kappa * (nth+1) + 12 * np.exp(-2*t**3) * (-np.sin(15*t)**2)"
    # gamma2 becomes negative during some time intervals

    # nm_mcsolve integration
    ops_and_rates = []
    ops_and_rates.append([a0.dag(), gamma1])
    ops_and_rates.append([a0,       gamma2])
    nm_options = {'map': 'parallel', 'improved_sampling': True}
    MCSol = nm_mcsolve(H, psi0, times, ops_and_rates,
                       args={'kappa': 1.0 / 0.129, 'nth': 0.063},
                       e_ops=[a0.dag() * a0, a0 * a0.dag()],
                       options=nm_options, ntraj=2500)

    # mesolve integration for comparison
    d_ops = [[lindblad_dissipator(a0.dag(), a0.dag()), gamma1],
             [lindblad_dissipator(a0, a0),             gamma2]]
    MESol = mesolve(H, psi0, times, d_ops, e_ops=[a0.dag() * a0, a0 * a0.dag()],
                    args={'kappa': 1.0 / 0.129, 'nth': 0.063})

    plt.figure()
    plt.plot(times, MCSol.expect[0], 'g',
             times, MCSol.expect[1], 'b',
             times, MCSol.trace, 'r')
    plt.plot(times, MESol.expect[0], 'g--',
             times, MESol.expect[1], 'b--')
    plt.title('Monte Carlo time evolution')
    plt.xlabel('Time')
    plt.ylabel('Expectation values')
    plt.legend((r'$\langle 1 | \rho | 1 \rangle$',
                r'$\langle 0 | \rho | 0 \rangle$',
                r'$\operatorname{tr} \rho$'))
    plt.show()