# Lecture 7 - Symbolic quantum mechanics using SymPsi - Semiclassical equations of motion¶

Author: J. R. Johansson ([email protected]), http://jrjohansson.github.io, and Eunjong Kim.

Status: Preliminary (work in progress)

This notebook is part of a series of IPython notebooks on symbolic quantum mechanics computations using SymPy and SymPsi. SymPsi is an experimental fork and extension of the sympy.physics.quantum module in SymPy. The latest version of this notebook is available at http://github.com/jrjohansson/sympy-quantum-notebooks, and the other notebooks in this lecture series are also indexed at http://jrjohansson.github.io.

Requirements: A recent version of SymPy and the latest development version of SymPsi is required to execute this notebook. Instructions for how to install SymPsi is available here.

Disclaimer: The SymPsi module is still under active development and may change in behavior without notice, and the intention is to move some of its features to sympy.physics.quantum when they matured and have been tested. However, these notebooks will be kept up-to-date the latest versions of SymPy and SymPsi.

## Setup modules¶

In [1]:
%matplotlib inline
import matplotlib.pyplot as plt

In [2]:
import numpy as np

In [3]:
from sympy import *
init_printing()

In [4]:
from sympsi import *
from sympsi.boson import *
from sympsi.pauli import *
from sympsi.operatorordering import *
from sympsi.expectation import *
from sympsi.operator import OperatorFunction


## Semiclassical equations of motion¶

The dynamics of an open quantum system with a given Hamiltonian, $H$, and some interaction with an environment that acts on the system through the sytem operator $a$, and with rate $\kappa$, can often be described with a Lindblad master equation for the dynamics of the system density matrix $\rho$:

$$\frac{d}{dt}\rho = -i[H, \rho] + \kappa \mathcal{D}[a]\rho,$$

where the Lindblad superoperator $\mathcal{D}$ is

$$\mathcal{D}[a]\rho = a \rho a^\dagger -\frac{1}{2}\rho a^\dagger a - \frac{1}{2}a^\dagger a \rho.$$

One common approach to solve for the dynamics of this system is to represent the system operators and the density operator as matrices, possibly in a truncated state space, and solve the matrix-valued ODE problem numerically.

Another approach is to use the adjoint master equation for the system operators $X$:

$$\frac{d}{dt} X = i [H, X] + \kappa \mathcal{D}[a^\dagger]X$$

and then solve for dynamics of the expectation values of the relevant system operators. The advantage of this method is that the ODEs are no longer matrix-valued, unlike the ODE for the density matrix. However, from the density matrix we can calculate any same-time expectation values, but with explicit ODEs for expectation values we need to select in advance which operator's expectation values we want to generate equations for.

We can easily generate an equation for the expectation value of a specific operator by multiplying the master equation for $\rho$ from the left with an operator $X$, and then take the trace over the entire equation. Doing this we obtain:

$$X\frac{d}{dt}\rho = -iX[H, \rho] + \kappa X\mathcal{D}[a]\rho$$

and taking the trace:

$${\rm Tr}\left(X\frac{d}{dt}\rho\right) = -i{\rm Tr}\left(X[H, \rho]\right) + \kappa {\rm Tr}\left(X\mathcal{D}[a]\rho\right)$$

using the cyclic permutation properties of traces:

$$\frac{d}{dt}{\rm Tr}\left(X\rho\right) = -i{\rm Tr}\left([X, H]\rho\right) + \kappa {\rm Tr}\left((\mathcal{D}[a]X) \rho\right)$$

we end up with an equation for the expectation value of the operator $X$:

$$# \frac{d}{dt}\langle X\rangle ¶ i\langle [H, X] \rangle + \kappa \langle \mathcal{D}[a]X \rangle$$

Note that this is a C-number equation, and therefore not as complicated to solve as the master equation for the density matrix. However, the problem with this C-number equation is that the expressions $[H, X]$ and $\mathcal{D}[a]X$ in general will introduce dependencies on other system operators, so we obtain a system of coupled C-number equations. If this system of equations closes when a finite number of operators are included, then we can use this method to solve the dynamics of these expectation values exactly. If the system of equations do not close, which is often the case for coupled systems, then we can still use this method if we introduce some rule for truncating high-order operator expectation values (for example, by discarding high-order terms or by factoring them in expectation values of lower order). However, in this case the results are no longer exact, and is called a semi-classical equation of motion.

With SymPsi we can automatically generate semiclassical equations of motion for operators in a system described by a given Hamiltonian and a set of collapse operators that describe its coupling to an environment.

## Driven harmonic oscillator¶

Consider a driven harmonic oscillator, which interaction with an bath at some temperature that corresponds to $N_{\rm th}$ average photons. We begin by setting up symbolic variables for the problem parameters and the system operators in SymPsi:

In [5]:
w, t, Nth, Ad, kappa = symbols(r"\omega, t, n_{th}, A_d, kappa", positive=True)

In [6]:
a = BosonOp("a")
rho = Operator(r"\rho")
rho_t = OperatorFunction(rho, t)

In [7]:
H = w * Dagger(a) * a + Ad * (a + Dagger(a))

Eq(Symbol("H"), H)

Out[7]:
$$H = A_{d} \left({{a}^\dagger} + {a}\right) + \omega {{a}^\dagger} {a}$$

The master equation for this system can be generated using the master_equation function:

In [8]:
c_ops = [sqrt(kappa * (Nth + 1)) * a, sqrt(kappa * Nth) * Dagger(a)]

In [9]:
me = master_equation(rho_t, t, H, c_ops)

me

Out[9]:
$$\frac{\partial}{\partial t} {{\rho}(t)} = \kappa n_{{th}} {{a}^\dagger} {{\rho}(t)} {a} - \frac{\kappa n_{{th}}}{2} {a} {{a}^\dagger} {{\rho}(t)} - \frac{\kappa n_{{th}}}{2} {{\rho}(t)} {a} {{a}^\dagger} - \frac{\kappa {{a}^\dagger}}{2} \left(n_{{th}} + 1\right) {a} {{\rho}(t)} + \kappa \left(n_{{th}} + 1\right) {a} {{\rho}(t)} {{a}^\dagger} - \frac{\kappa {{\rho}(t)}}{2} \left(n_{{th}} + 1\right) {{a}^\dagger} {a} - i \left[A_{d} \left({{a}^\dagger} + {a}\right) + \omega {{a}^\dagger} {a},{{\rho}(t)}\right]$$

Equation for the system operators can be generated using the function operator_master_equation, and for the specific case of the cavity operator $a$ we obtain:

In [10]:
# first setup time-dependent operators
a_t = OperatorFunction(a, t)
a_to_a_t = {a: a_t, Dagger(a): Dagger(a_t)}
H_t = H.subs(a_to_a_t)
c_ops_t = [c.subs(a_to_a_t) for c in c_ops]

In [11]:
# operator master equation for a
ome_a = operator_master_equation(a_t, t, H_t, c_ops_t)

Eq(ome_a.lhs, normal_ordered_form(ome_a.rhs.doit().expand()))

Out[11]:
$$\frac{d}{d t} {{{a}}(t)} = - i A_{d} - i \omega {{{a}}(t)} - \frac{\kappa {{{a}}(t)}}{2}$$
In [12]:
# operator master equation for n = Dagger(a) * a
ome_n = operator_master_equation(Dagger(a_t) * a_t, t, H_t, c_ops_t)

Eq(ome_n.lhs, normal_ordered_form(ome_n.rhs.doit().expand()))

Out[12]:
$$\frac{d}{d t} {{{{a}^\dagger}}(t)} {{{a}}(t)} + {{{{a}^\dagger}}(t)} \frac{d}{d t} {{{a}}(t)} = - i A_{d} {{{{a}^\dagger}}(t)} + i A_{d} {{{a}}(t)} + \kappa n_{{th}} - \kappa {{{{a}^\dagger}}(t)} {{{a}}(t)}$$

From these operator equations we see that the equation for $a$ depends only on the operator $a$, while the equation for $n$ depends on $n$, $a$, $a^\dagger$. So to solve the latter equation we therefore also have to generate an equation for $a^\dagger$.

### System of semiclassical equations¶

In [13]:
ops, op_eqm, sc_eqm, sc_ode, ofm, oim = semi_classical_eqm(H, c_ops)

In [14]:
html_table([[Eq(Expectation(key), ofm[key]), sc_ode[key]] for key in operator_sort_by_order(sc_ode)])

Out[14]:
 $\left\langle {{a}^\dagger} \right\rangle = \operatorname{A_{0}}{\left (t \right )}$ $\frac{d}{d t} \operatorname{A_{0}}{\left (t \right )} = i A_{d} + i \omega \operatorname{A_{0}}{\left (t \right )} - \frac{\kappa}{2} \operatorname{A_{0}}{\left (t \right )}$ $\left\langle {a} \right\rangle = \operatorname{A_{1}}{\left (t \right )}$ $\frac{d}{d t} \operatorname{A_{1}}{\left (t \right )} = - i A_{d} - i \omega \operatorname{A_{1}}{\left (t \right )} - \frac{\kappa}{2} \operatorname{A_{1}}{\left (t \right )}$ $\left\langle {{a}^\dagger} {a} \right\rangle = \operatorname{A_{2}}{\left (t \right )}$ $\frac{d}{d t} \operatorname{A_{2}}{\left (t \right )} = - i A_{d} \operatorname{A_{0}}{\left (t \right )} + i A_{d} \operatorname{A_{1}}{\left (t \right )} + \kappa n_{{th}} - \kappa \operatorname{A_{2}}{\left (t \right )}$

Since this is a system of linear ODEs, we can write it on matrix form:

In [15]:
A_eq, A, M, b = semi_classical_eqm_matrix_form(sc_ode, t, ofm)

A_eq

Out[15]:
$$- \frac{d}{d t} \left[\begin{matrix}\operatorname{A_{0}}{\left (t \right )}\\\operatorname{A_{1}}{\left (t \right )}\\\operatorname{A_{2}}{\left (t \right )}\end{matrix}\right] = \left[\begin{matrix}i A_{d}\\- i A_{d}\\\kappa n_{{th}}\end{matrix}\right] + \left[\begin{matrix}i \omega - \frac{\kappa}{2} & 0 & 0\\0 & - i \omega - \frac{\kappa}{2} & 0\\- i A_{d} & i A_{d} & - \kappa\end{matrix}\right] \left[\begin{matrix}\operatorname{A_{0}}{\left (t \right )}\\\operatorname{A_{1}}{\left (t \right )}\\\operatorname{A_{2}}{\left (t \right )}\end{matrix}\right]$$

We can solve for the steadystate by setting the left-hand-side of the ODE to zero, and solve the linear system of equations:

In [16]:
A_sol = M.LUsolve(-b)


The solution for the three system operators are:

In [17]:
A_sol[ops.index(Dagger(a)*a)]

Out[17]:
$$- \frac{1}{\kappa} \left(\frac{A_{d}^{2}}{i \omega - \frac{\kappa}{2}} + \frac{A_{d}^{2}}{- i \omega - \frac{\kappa}{2}} - \kappa n_{{th}}\right)$$
In [18]:
A_sol[ops.index(a)]

Out[18]:
$$\frac{i A_{d}}{- i \omega - \frac{\kappa}{2}}$$
In [19]:
A_sol[ops.index(Dagger(a))]

Out[19]:
$$- \frac{i A_{d}}{i \omega - \frac{\kappa}{2}}$$

We can also solve for the steadystate directly from the ODE by settings its right-hand-side to zero, and using the SymPy solve function:

In [20]:
solve([eq.rhs for eq in sc_ode.values()], list(ofm.values()))

Out[20]:
$$\left \{ \operatorname{A_{0}}{\left (t \right )} : - \frac{2 i A_{d}}{2 i \omega - \kappa}, \quad \operatorname{A_{1}}{\left (t \right )} : - \frac{2 i A_{d}}{2 i \omega + \kappa}, \quad \operatorname{A_{2}}{\left (t \right )} : \frac{1}{4 \omega^{2} + \kappa^{2}} \left(4 A_{d}^{2} + 4 \omega^{2} n_{{th}} + \kappa^{2} n_{{th}}\right)\right \}$$

### Solve in the ODEs¶

For systems with a small number of dependent operators we can solve the resulting system of ODEs directly:

In [21]:
sols = dsolve(list(sc_ode.values())); sols

Out[21]:
$$\left [ \operatorname{A_{0}}{\left (t \right )} = C_{1} e^{t \left(i \omega - \frac{\kappa}{2}\right)}, \quad \operatorname{A_{1}}{\left (t \right )} = C_{2} e^{t \left(- i \omega - \frac{\kappa}{2}\right)}, \quad \operatorname{A_{2}}{\left (t \right )} = - \frac{i A_{d} C_{1} e^{t \left(i \omega - \frac{\kappa}{2}\right)}}{i \omega + \frac{\kappa}{2}} + \frac{i A_{d} C_{2} e^{t \left(- i \omega - \frac{\kappa}{2}\right)}}{- i \omega + \frac{\kappa}{2}} + \frac{C_{3}}{e^{\kappa t}}\right ]$$
In [22]:
# hack
tt = [s for s in sols[0].rhs.free_symbols if s.name == 't'][0]


We also need o specify the initial conditions: Here the initial conditions are $\langle a(0) \rangle = \langle a^\dagger(0) \rangle = 2$ and $\langle a^\dagger(0)a(0) \rangle = 4$.

In [23]:
ics = {ofm[Dagger(a)].subs(tt, 0): 2,
ofm[a].subs(tt, 0): 2,
ofm[Dagger(a)*a].subs(tt, 0): 4}; ics

Out[23]:
$$\left \{ \operatorname{A_{0}}{\left (0 \right )} : 2, \quad \operatorname{A_{1}}{\left (0 \right )} : 2, \quad \operatorname{A_{2}}{\left (0 \right )} : 4\right \}$$
In [24]:
constants = set(sum([[s for s in sol.free_symbols if (str(s)[0] == 'C')] for sol in sols], [])); constants

Out[24]:
$$\left\{C_{1}, C_{2}, C_{3}\right\}$$
In [25]:
C_sols = solve([sol.subs(tt, 0).subs(ics) for sol in sols], constants); C_sols

Out[25]:
$$\left \{ C_{1} : 2, \quad C_{2} : 2, \quad C_{3} : \frac{1}{4 \omega^{2} + \kappa^{2}} \left(16 A_{d} \omega + 16 \omega^{2} + 4 \kappa^{2}\right)\right \}$$
In [26]:
sols_with_ics = [sol.subs(C_sols) for sol in sols]; sols_with_ics

Out[26]:
$$\left [ \operatorname{A_{0}}{\left (t \right )} = 2 e^{t \left(i \omega - \frac{\kappa}{2}\right)}, \quad \operatorname{A_{1}}{\left (t \right )} = 2 e^{t \left(- i \omega - \frac{\kappa}{2}\right)}, \quad \operatorname{A_{2}}{\left (t \right )} = - \frac{2 i A_{d} e^{t \left(i \omega - \frac{\kappa}{2}\right)}}{i \omega + \frac{\kappa}{2}} + \frac{2 i A_{d} e^{t \left(- i \omega - \frac{\kappa}{2}\right)}}{- i \omega + \frac{\kappa}{2}} + \frac{16 A_{d} \omega + 16 \omega^{2} + 4 \kappa^{2}}{\left(4 \omega^{2} + \kappa^{2}\right) e^{\kappa t}}\right ]$$

Now let's insert numerical values for the system parameters so we can plot the solution:

In [27]:
values = {w: 1.0, Ad: 0.0, kappa: 0.1, Nth: 0.0}

In [28]:
sols_funcs = [sol.rhs.subs(values) for sol in sols_with_ics]; sols_funcs

Out[28]:
$$\left [ 2 e^{t \left(-0.05 + 1.0 i\right)}, \quad 2 e^{t \left(-0.05 - 1.0 i\right)}, \quad \frac{4}{e^{0.1 t}}\right ]$$
In [29]:
times = np.linspace(0, 50, 500)

y_funcs = [lambdify([tt], sol_func, 'numpy') for sol_func in sols_funcs]

fig, axes = plt.subplots(len(y_funcs), 1, figsize=(12, 6))

for n, y_func in enumerate(y_funcs):
axes[n].plot(times, np.real(y_func(times)), 'r')

axes[2].set_ylim(0, 5);


## Driven dissipative two-level system¶

In [30]:
sx, sy, sz, sm, sp = SigmaX(), SigmaY(), SigmaZ(), SigmaMinus(), SigmaPlus()

In [31]:
Omega, gamma_0, N, t = symbols("\Omega, \gamma_0, N, t", positive=True)

values = {Omega: 1.0, gamma_0: 0.5, N: 1.75}

In [32]:
H = -Omega/2 * sx
H

Out[32]:
$$- \frac{\Omega {\sigma_x}}{2}$$
In [33]:
c_ops = [sqrt(gamma_0 * (N + 1)) * pauli_represent_x_y(sm),
sqrt(gamma_0 * N) * pauli_represent_x_y(sp)]

In [34]:
ops, op_eqm, sc_eqm, sc_ode, ofm, oim = semi_classical_eqm(H, c_ops)

In [35]:
html_table([[Eq(Expectation(key), ofm[key]), sc_ode[key]] for key in operator_sort_by_order(sc_ode)])

Out[35]:
 $\left\langle {\sigma_x} \right\rangle = \operatorname{A_{0}}{\left (t \right )}$ $\frac{d}{d t} \operatorname{A_{0}}{\left (t \right )} = - N \gamma_{0} \operatorname{A_{0}}{\left (t \right )} - \frac{\gamma_{0}}{2} \operatorname{A_{0}}{\left (t \right )}$ $\left\langle {\sigma_y} \right\rangle = \operatorname{A_{1}}{\left (t \right )}$ $\frac{d}{d t} \operatorname{A_{1}}{\left (t \right )} = - N \gamma_{0} \operatorname{A_{1}}{\left (t \right )} + \Omega \operatorname{A_{2}}{\left (t \right )} - \frac{\gamma_{0}}{2} \operatorname{A_{1}}{\left (t \right )}$ $\left\langle {\sigma_z} \right\rangle = \operatorname{A_{2}}{\left (t \right )}$ $\frac{d}{d t} \operatorname{A_{2}}{\left (t \right )} = - 2 N \gamma_{0} \operatorname{A_{2}}{\left (t \right )} - \Omega \operatorname{A_{1}}{\left (t \right )} - \gamma_{0} \operatorname{A_{2}}{\left (t \right )} - \gamma_{0}$
In [36]:
A_eq, A, M, b = semi_classical_eqm_matrix_form(sc_ode, t, ofm)

A_eq

Out[36]:
$$- \frac{d}{d t} \left[\begin{matrix}\operatorname{A_{0}}{\left (t \right )}\\\operatorname{A_{1}}{\left (t \right )}\\\operatorname{A_{2}}{\left (t \right )}\end{matrix}\right] = \left[\begin{matrix}0\\0\\- \gamma_{0}\end{matrix}\right] + \left[\begin{matrix}- N \gamma_{0} - \frac{\gamma_{0}}{2} & 0 & 0\\0 & - N \gamma_{0} - \frac{\gamma_{0}}{2} & \Omega\\0 & - \Omega & - 2 N \gamma_{0} - \gamma_{0}\end{matrix}\right] \left[\begin{matrix}\operatorname{A_{0}}{\left (t \right )}\\\operatorname{A_{1}}{\left (t \right )}\\\operatorname{A_{2}}{\left (t \right )}\end{matrix}\right]$$

In [37]:
A_sol = M.LUsolve(-b)


The steadystate expectation value of $\sigma_x$:

In [38]:
A_sol[ops.index(sx)]

Out[38]:
$$0$$

The steadystate expectation value of $\sigma_y$:

In [39]:
A_sol[ops.index(sy)]

Out[39]:
$$- \frac{\Omega \gamma_{0}}{\left(- N \gamma_{0} - \frac{\gamma_{0}}{2}\right) \left(- 2 N \gamma_{0} + \frac{\Omega^{2}}{- N \gamma_{0} - \frac{\gamma_{0}}{2}} - \gamma_{0}\right)}$$

The steadystate expectation value of $\sigma_z$:

In [40]:
A_sol[ops.index(sz)]

Out[40]:
$$\frac{\gamma_{0}}{- 2 N \gamma_{0} + \frac{\Omega^{2}}{- N \gamma_{0} - \frac{\gamma_{0}}{2}} - \gamma_{0}}$$

Steadystate of $\sigma_+$:

In [41]:
pauli_represent_x_y(sp).subs({sx: A_sol[ops.index(sx)], sy: A_sol[ops.index(sy)]})

Out[41]:
$$- \frac{i \Omega \gamma_{0}}{2 \left(- N \gamma_{0} - \frac{\gamma_{0}}{2}\right) \left(- 2 N \gamma_{0} + \frac{\Omega^{2}}{- N \gamma_{0} - \frac{\gamma_{0}}{2}} - \gamma_{0}\right)}$$

Steadystate of $\sigma_-$:

In [42]:
pauli_represent_x_y(sm).subs({sx: A_sol[ops.index(sx)], sy: A_sol[ops.index(sy)]})

Out[42]:
$$\frac{i \Omega \gamma_{0}}{2 \left(- N \gamma_{0} - \frac{\gamma_{0}}{2}\right) \left(- 2 N \gamma_{0} + \frac{\Omega^{2}}{- N \gamma_{0} - \frac{\gamma_{0}}{2}} - \gamma_{0}\right)}$$

Alternatively we can also use the SymPy solve function to find the steadystate solutions:

In [43]:
solve([eq.rhs for eq in sc_ode.values()], list(ofm.values()))

Out[43]:
$$\left \{ \operatorname{A_{0}}{\left (t \right )} : 0, \quad \operatorname{A_{1}}{\left (t \right )} : - \frac{2 \Omega \gamma_{0}}{2 \Omega^{2} + \gamma_{0}^{2} \left(2 N + 1\right)^{2}}, \quad \operatorname{A_{2}}{\left (t \right )} : - \frac{\gamma_{0}^{2} \left(2 N + 1\right)}{2 \Omega^{2} + \gamma_{0}^{2} \left(2 N + 1\right)^{2}}\right \}$$

At zero temperature:

In [44]:
solve([eq.subs(N, 0).rhs for eq in sc_ode.values()], list(ofm.values()))

Out[44]:
$$\left \{ \operatorname{A_{0}}{\left (t \right )} : 0, \quad \operatorname{A_{1}}{\left (t \right )} : - \frac{2 \Omega \gamma_{0}}{2 \Omega^{2} + \gamma_{0}^{2}}, \quad \operatorname{A_{2}}{\left (t \right )} : - \frac{\gamma_{0}^{2}}{2 \Omega^{2} + \gamma_{0}^{2}}\right \}$$

## Versions¶

In [45]:
%reload_ext version_information

%version_information sympy, sympsi

Out[45]:
SoftwareVersion
Python3.4.1 (default, Sep 20 2014, 19:44:17) [GCC 4.2.1 Compatible Apple LLVM 5.1 (clang-503.0.40)]
IPython2.3.0
OSDarwin 13.4.0 x86_64 i386 64bit
sympy0.7.5-git
sympsi0.1.0.dev-0c6e514
Sun Oct 12 21:42:08 2014 JST