Reference

Question

In the book Breuer, Heinz-Peter, and Francesco Petruccione. The Theory of Open Quantum Systems. Oxford: Clarendon Press, 2009., there is equation(10.17)

\[ \frac{\mathrm{d}}{\mathrm{d}t}c_1(t) = - \int_0^t \mathrm{d}t_1 f(t-t_1) c_1(t_1) \]

where

\[ f(t) = \frac{1}{2} \gamma_0\lambda e^{-\lambda |t|} \]

\(\gamma_0\) and \(\lambda\) are parameters, for example, we can set \(\gamma_0=1, \lambda=10\). For this \(f(t)\), we can solve it analytically, the solution is

\[ c_1(t) = c_1(0) e^{-\lambda t/2} \left[\cosh\left(\frac{dt}{2}\right) + \frac{\lambda}{d} \sinh\left(\frac{dt}{2}\right) \right] \]

where \(d = \sqrt{\lambda^2 - 2\gamma_0\lambda}\).

My question is, how to numerically solve (for example, using python function scipy.integrate.solve_ivp, which is a function solving OEDs by Runge-Kutta method or other methods) this equation?

I want to numerically solve it, because

In [1]

import numpy as np
from scipy.integrate import quad, solve_ivp
import qutip as qutip
import matplotlib.pyplot as plt

In [2]

def ft(t, g0=1, lam=10):
    """Breuer Eq.(10.44)"""
    return 1/2 * g0 * lam * np.exp(-lam*np.abs(t))


def dc_dt(t, c0=1, lam=10, g0=1):
    """Right hand side of Breuer Eq.(10.17),
    chose f(t-t1) as  Breuer Eq.(10.44)"""
    res = quad(lambda t1: ft(t-t1)*ct(t1), 0, t)[0]
    return -res


def ct(t, c0=1, lam=10, g0=1):
    """Breuer Eq.(10.45)"""
    d = np.emath.sqrt(lam**2 - 2*g0*lam)
    ct = c0 * np.exp(-lam*t/2)
    ct *= np.cosh(d*t/2) + lam/d*np.sinh(d*t/2)
    return ct


def ct0(t, c0=1, g0=1):
    """Markovian limit solution"""
    return c0*np.exp(-g0*t/2)

In [3]

t = np.linspace(0, 10, 1000)
plt.plot(t, [np.abs(ct(ti, lam=1)) for ti in t], 'y')
plt.plot(t, np.abs(ct(t, lam=2.1)), 'r')
plt.plot(t, np.abs(ct(t, lam=10)), 'g')
plt.plot(t, np.abs(ct(t, lam=100)), 'b')
plt.plot(t[::10], ct0(t)[::10], 'kx', label='Markov')
plt.xlabel(r'$t$')
plt.legend()
<Figure size 640x480 with 1 Axes>

In [4]

# check Breuer Eq.(10.45)
dt = 1e-2
t = np.arange(0, 4, dt)


def ft(t, g0=1, lam=10):
    """Breuer Eq.(10.44)"""
    return 1/2 * g0 * lam * np.exp(-lam*np.abs(t))


def dc_dt(t, c0=1, lam=10, g0=1):
    """Right hand side of Breuer Eq.(10.17),
    chose f(t-t1) as  Breuer Eq.(10.44)"""
    res = quad(lambda t1: ft(t-t1)*ct(t1), 0, t)[0]
    return -res

cts = ct(t)
plt.plot(t[1:], np.diff(cts)/dt)

plt.plot(t, np.array([dc_dt(ti) for ti in t]), 'x')
<Figure size 640x480 with 1 Axes>

In [5]

dt = 1e-2
ts = np.arange(0, 5, dt)
dcdt = [dc_dt(ti) for ti in ts]
plt.plot(ts, dcdt, 'r--')
plt.plot(ts[1:], np.diff(dcdt)/dt, 'g')
def myd(t):
    return - ft(0)*ct(t) + ft(t)*ct(0)
plt.plot(ts, myd(ts), 'k')


<Figure size 640x480 with 1 Axes>

In [2]

print('hello')
hello