# Lecture 3 - Symbolic quantum mechanics using SymPsi - Resonators and cavities¶

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

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]:
from sympy import *
init_printing()

In [2]:
from sympsi import *
from sympsi.boson import *
from sympsi.operatorordering import *


## Introduction¶

In this notebook we will work with cavities and resonators. A single mode of a resonator can be modelled with a quantum harmonic oscillator. Here we look at the effect of classical driving fields, transformation to different rotating frames, and various types of coupling between different modes in a resonator.

## Quantum Harmonic Oscillator¶

First we consider the Hamiltonian for a single harmonic oscillator, which describes a single mode of for example a cavity or a waveguide resonator:

$$H = \hbar \omega_r a^\dagger a$$

To represent this Hamiltonian in sympy we create a symbol omega_r and an instance of the class BosonOp:

In [3]:
omega_r = symbols("omega_r", positive=True)
Hsym = symbols("H")
a = BosonOp("a")

In [4]:
H0 = omega_r * Dagger(a) * a

Eq(Hsym, H0)

Out[4]:
$$H = \omega_{r} {{a}^\dagger} {a}$$

### Classical drive signal¶

A classical driving field of frequency $\omega_d$ and amplitude $A$ and phase $\phi_0$ can be modelled by including an additional term in the Hamiltonian:

$$H = \hbar \omega_r a^\dagger a + A \cos(\omega_dt + \phi_0) (a + a^\dagger)$$

where we have assumed that the driving field couples to the quadrature of the resonator $a + a^\dagger$, which is the canonical case. It is convenient to rewrite the $\cos$ factor

$$H = \hbar \omega_r a^\dagger a + (Ae^{-i\omega_dt}+ A^*e^{i\omega_dt})(a + a^\dagger)$$

where we have redefined $A \rightarrow Ae^{-i\phi_0}$. In Sympy we can represent this Hamiltonian as:

In [5]:
omega_d, t = symbols("omega_d, t")
A = symbols("A")

In [6]:
Hdrive = (A * exp(-I * omega_d * t) + conjugate(A) * exp(I * omega_d * t)) * (a + Dagger(a))

Hdrive

Out[6]:
$$\left(A e^{- i \omega_{d} t} + e^{i \omega_{d} t} \overline{A}\right) \left({{a}^\dagger} + {a}\right)$$
In [7]:
H = H0 + Hdrive

Eq(Hsym, H)

Out[7]:
$$H = \omega_{r} {{a}^\dagger} {a} + \left(A e^{- i \omega_{d} t} + e^{i \omega_{d} t} \overline{A}\right) \left({{a}^\dagger} + {a}\right)$$

When working with Hamiltonians like this one, one common operation is to move to different rotating frames (by performing unitary transformations) where the Hamiltonian takes a simplier form.

Here we want to transform this Hamiltonian to a rotating frame in which the drive term (and all other terms) are no longer explicitly depening on time t. We can accomplish this by performing the unitary transformation

$$U = \exp(i\omega_d a^\dagger a t)$$

In [8]:
U = exp(I * omega_d * t * Dagger(a) * a)

U

Out[8]:
$$e^{i \omega_{d} t {{a}^\dagger} {a}}$$

When doing a unitary basis transformation, the Hamiltonian is transformed as

$$H \rightarrow UHU^\dagger -i U \frac{d}{dt}U^\dagger$$

and we can carry out this transformation using the function hamiltonian_transformation:

In [9]:
H2 = hamiltonian_transformation(U, H.expand())

H2

Out[9]:
$$A {{a}^\dagger} + A e^{- 2 i \omega_{d} t} {a} - \omega_{d} {{a}^\dagger} {a} + \omega_{r} {{a}^\dagger} {a} + e^{2 i \omega_{d} t} \overline{A} {{a}^\dagger} + \overline{A} {a}$$

Now let's introduce the detuning $\Delta = \omega_r - \omega_d$:

In [10]:
Delta = symbols("Delta", positive=True)

In [11]:
H3 = collect(H2, Dagger(a) * a).subs(omega_r - omega_d, Delta)

H3

Out[11]:
$$A {{a}^\dagger} + A e^{- 2 i \omega_{d} t} {a} + \Delta {{a}^\dagger} {a} + e^{2 i \omega_{d} t} \overline{A} {{a}^\dagger} + \overline{A} {a}$$

### Rotating-wave approximation (RWA)¶

Now we invoke the rotating-wave approximation, under which we assume that $\omega_d$ is much larger than $\Delta$, so that the we can neglect the two fast rotating terms which contains factors $\exp(\pm 2i \omega_d t)$:

In [12]:
H3 = drop_terms_containing(H3, [exp( 2 * I * omega_d * t),
exp(-2 * I * omega_d * t)])

Eq(Hsym, H3)

Out[12]:
$$H = A {{a}^\dagger} + \Delta {{a}^\dagger} {a} + \overline{A} {a}$$

This is a time-independent hamiltonian describing the resonator in a rotating frame, where fast rotating terms have been dropped. We can no apply a displacement transformation by applying the unitary transformation:

In [13]:
alpha = symbols("alpha")
H = Dagger(a) * alpha - conjugate(alpha) * a
U = exp(H)

U

Out[13]:
$$e^{\alpha {{a}^\dagger} - \overline{\alpha} {a}}$$
In [14]:
H4 = hamiltonian_transformation(U, H3)

H4 = collect(H4.expand(), [Dagger(a)*a, a, Dagger(a)])

H4

Out[14]:
$$- A \overline{\alpha} + \Delta \alpha \overline{\alpha} + \Delta {{a}^\dagger} {a} - \alpha \overline{A} + \left(A - \Delta \alpha\right) {{a}^\dagger} + \left(- \Delta \overline{\alpha} + \overline{A}\right) {a}$$

If we choose $\alpha = A/\Delta$ and drop c-numbers from the hamiltonian, we obtain a particularly simple form:

In [15]:
H5 = H4.subs(alpha, A/Delta)

H5 = drop_c_number_terms(H5)

H5

Out[15]:
$$\Delta {{a}^\dagger} {a}$$

So a driven harmonic oscillator can be described by the hamiltonian of an undriven harmonic oscillator in a displaced frame.

## Optical parametric oscillator¶

An optical parametric oscillator (OPO) is a two-mode system with a particular nonlinear interaction. The Hamiltonian of the OPO is

$$H = \omega_r a^\dagger a + \omega_p b^\dagger b + (\kappa {a^\dagger}^2 b + \kappa^* a^2 b^\dagger) + (A e^{i\omega_pt} + A^* e^{-i\omega_pt})(b + b^\dagger)$$

where the operators of the two modes are$a$ and $b$, and the nonlinear interaction strength is $\kappa$, and where we have also included a classical drive field applied to mode $b$.

To model this system in SymPy we create two instances of BosonOp, one for each mode, and construct the Hamiltonian:

In [16]:
a, b = BosonOp("a"), BosonOp("b")

In [17]:
kappa, omega_p, omega_a, omega_b, t = symbols("kappa, omega_p, omega_a, omega_b, t")
Delta_a, Delta_b = symbols("Delta_a, Delta_b", positive=True)

In [18]:
Hdrive = (A * exp(-I * omega_p * t) + conjugate(A) * exp(I * omega_p * t)) * (b + Dagger(b))
H1 = omega_a * Dagger(a) * a + omega_b * Dagger(b) * b  + kappa * (a ** 2 * Dagger(b) + Dagger(a) ** 2 * b) + Hdrive

Eq(Hsym, H1)

Out[18]:
$$H = \kappa \left(\left({{a}^\dagger}\right)^{2} {b} + \left({a}\right)^{2} {{b}^\dagger}\right) + \omega_{a} {{a}^\dagger} {a} + \omega_{b} {{b}^\dagger} {b} + \left(A e^{- i \omega_{p} t} + e^{i \omega_{p} t} \overline{A}\right) \left({{b}^\dagger} + {b}\right)$$

We first move to a frame rotating with frequency $\omega_p$ with respect to mode $b$:

In [19]:
U = exp(I * omega_p * t * Dagger(b) * b)

U

Out[19]:
$$e^{i \omega_{p} t {{b}^\dagger} {b}}$$
In [20]:
H2 = hamiltonian_transformation(U, H1.expand(), independent=True)

H2

Out[20]:
$$A {{b}^\dagger} + A e^{- 2 i \omega_{p} t} {b} + \kappa e^{i \omega_{p} t} \left({a}\right)^{2} {{b}^\dagger} + \kappa e^{- i \omega_{p} t} \left({{a}^\dagger}\right)^{2} {b} + \omega_{a} {{a}^\dagger} {a} + \omega_{b} {{b}^\dagger} {b} - \omega_{p} {{b}^\dagger} {b} + e^{2 i \omega_{p} t} \overline{A} {{b}^\dagger} + \overline{A} {b}$$

and we can perform the rotating wave approximation and drop the fast rotating terms $\exp{\left(\pm i \omega_p t\right)}$:

In [21]:
H3 = drop_terms_containing(H2, [exp(2 * I * omega_p * t),
exp(-2 * I * omega_p * t)])

H3 = collect(H3, kappa)

H3

Out[21]:
$$A {{b}^\dagger} + \kappa \left(e^{i \omega_{p} t} \left({a}\right)^{2} {{b}^\dagger} + e^{- i \omega_{p} t} \left({{a}^\dagger}\right)^{2} {b}\right) + \omega_{a} {{a}^\dagger} {a} + \omega_{b} {{b}^\dagger} {b} - \omega_{p} {{b}^\dagger} {b} + \overline{A} {b}$$

Introduce a new variable for the detuning $\Delta_b = \omega_b - \omega_p$:

In [22]:
H4 = H3.subs(omega_b, Delta_b + omega_p).expand()

H4

Out[22]:
$$A {{b}^\dagger} + \Delta_{b} {{b}^\dagger} {b} + \kappa e^{i \omega_{p} t} \left({a}\right)^{2} {{b}^\dagger} + \kappa e^{- i \omega_{p} t} \left({{a}^\dagger}\right)^{2} {b} + \omega_{a} {{a}^\dagger} {a} + \overline{A} {b}$$

Next we want to displace the mode $b$ so that the drive terms are eliminated:

In [23]:
beta = symbols("beta")
H = Dagger(b) * beta - conjugate(beta) * b
U = exp(H)

U

Out[23]:
$$e^{\beta {{b}^\dagger} - \overline{\beta} {b}}$$
In [24]:
H5 = hamiltonian_transformation(U, H4.expand(), independent=True)

H5

Out[24]:
$$A \left(- \overline{\beta} + {{b}^\dagger}\right) + \Delta_{b} \left(- \overline{\beta} + {{b}^\dagger}\right) \left(- \beta + {b}\right) + \kappa e^{i \omega_{p} t} \left({a}\right)^{2} \left(- \overline{\beta} + {{b}^\dagger}\right) + \kappa e^{- i \omega_{p} t} \left({{a}^\dagger}\right)^{2} \left(- \beta + {b}\right) + \omega_{a} {{a}^\dagger} {a} + \overline{A} \left(- \beta + {b}\right)$$
In [25]:
H5 = collect(H5.expand(), [Dagger(a) * a, Dagger(b) * b, Dagger(a) ** 2 * b, a ** 2 * Dagger(b), b, Dagger(b)])

H5

Out[25]:
$$- A \overline{\beta} + \Delta_{b} \beta \overline{\beta} + \Delta_{b} {{b}^\dagger} {b} - \beta \kappa e^{- i \omega_{p} t} \left({{a}^\dagger}\right)^{2} - \beta \overline{A} - \kappa e^{i \omega_{p} t} \overline{\beta} \left({a}\right)^{2} + \kappa e^{i \omega_{p} t} \left({a}\right)^{2} {{b}^\dagger} + \kappa e^{- i \omega_{p} t} \left({{a}^\dagger}\right)^{2} {b} + \omega_{a} {{a}^\dagger} {a} + \left(A - \Delta_{b} \beta\right) {{b}^\dagger} + \left(- \Delta_{b} \overline{\beta} + \overline{A}\right) {b}$$

The choice $\beta = A/\Delta_b$ eliminates the drive terms. After dropping c numbers we have:

In [26]:
H6 = H5.subs(beta, A/Delta_b)

H6 = drop_c_number_terms(H6)

H6 = collect(H6, [exp( I * omega_p * t) * kappa * a ** 2,
exp(-I * omega_p * t) * kappa * Dagger(a) ** 2])

H6

Out[26]:
$$\Delta_{b} {{b}^\dagger} {b} + \kappa e^{i \omega_{p} t} \left({a}\right)^{2} \left({{b}^\dagger} - \frac{\overline{A}}{\Delta_{b}}\right) + \kappa e^{- i \omega_{p} t} \left({{a}^\dagger}\right)^{2} \left(- \frac{A}{\Delta_{b}} + {b}\right) + \omega_{a} {{a}^\dagger} {a}$$

If we now assume that the dynamics of the mode $b$ is dominated by the classical drive field, we can neglect the $b$ operators in the hamiltonian (in this displaced frame, which describes deviations from the dynamics induced by the classical driving field):

In [27]:
H7 = drop_terms_containing(H6.expand(), [b, Dagger(b)])

H7 = collect(H7, kappa / Delta_b)

H7

Out[27]:
$$\omega_{a} {{a}^\dagger} {a} + \frac{\kappa}{\Delta_{b}} \left(- A e^{- i \omega_{p} t} \left({{a}^\dagger}\right)^{2} - e^{i \omega_{p} t} \overline{A} \left({a}\right)^{2}\right)$$

To simpify the notation we redefine $- \kappa A / \Delta_b \rightarrow \kappa$, and assume that $A$ is real:

In [28]:
H = omega_a * Dagger(a) * a  + kappa * (a ** 2 * exp(I * omega_p * t) + Dagger(a) ** 2 * exp(-I * omega_p * t))

Eq(Hsym, H)

Out[28]:
$$H = \kappa \left(e^{i \omega_{p} t} \left({a}\right)^{2} + e^{- i \omega_{p} t} \left({{a}^\dagger}\right)^{2}\right) + \omega_{a} {{a}^\dagger} {a}$$

To eliminate the time-dependence in the interaction term we move to a rotating frame using the unitary transformation:

In [29]:
U = exp(I * omega_d * t * Dagger(a) * a)

U

Out[29]:
$$e^{i \omega_{d} t {{a}^\dagger} {a}}$$
In [30]:
H2 = hamiltonian_transformation(U, H)

H2

Out[30]:
$$\kappa \left(e^{2 i \omega_{d} t} e^{- i \omega_{p} t} \left({{a}^\dagger}\right)^{2} + e^{- 2 i \omega_{d} t} e^{i \omega_{p} t} \left({a}\right)^{2}\right) + \omega_{a} {{a}^\dagger} {a} - \omega_{d} {{a}^\dagger} {a}$$

Now consider the case $\omega_p = 2 \omega_d$:

In [31]:
H3 = H2.subs(omega_d, omega_p/2)

H3 = collect(H3, Dagger(a) * a)

H3

Out[31]:
$$\kappa \left(\left({{a}^\dagger}\right)^{2} + \left({a}\right)^{2}\right) + \left(\omega_{a} - \frac{\omega_{p}}{2}\right) {{a}^\dagger} {a}$$

Introduce $\Delta_a = \omega_a - \omega_p / 2$:

In [32]:
H4 = H3.subs(omega_a, Delta_a + omega_p/2)

H4

Out[32]:
$$\Delta_{a} {{a}^\dagger} {a} + \kappa \left(\left({{a}^\dagger}\right)^{2} + \left({a}\right)^{2}\right)$$

We can diagonalize this Hamiltonian by introducing the squeezing transformation:

In [33]:
chi = symbols("chi")

In [34]:
U = exp(chi/2 * a **2 - chi/2 * Dagger(a) ** 2)

U

Out[34]:
$$e^{- \frac{\chi \left({{a}^\dagger}\right)^{2}}{2} + \frac{\chi \left({a}\right)^{2}}{2}}$$

which tranforms the operators $a$ and $a^\dagger$ according to the well-known relations:

In [35]:
hamiltonian_transformation(U, a)

Out[35]:
$$\sinh{\left (\chi \right )} {{a}^\dagger} + \cosh{\left (\chi \right )} {a}$$
In [36]:
hamiltonian_transformation(U, Dagger(a))

Out[36]:
$$\sinh{\left (\chi \right )} {a} + \cosh{\left (\chi \right )} {{a}^\dagger}$$

and in this frame the Hamiltonian takes the form:

In [37]:
H5 = hamiltonian_transformation(U, H4)

H5

Out[37]:
$$\Delta_{a} \left(\sinh{\left (\chi \right )} {a} + \cosh{\left (\chi \right )} {{a}^\dagger}\right) \left(\sinh{\left (\chi \right )} {{a}^\dagger} + \cosh{\left (\chi \right )} {a}\right) + \kappa \left(\left(\sinh{\left (\chi \right )} {{a}^\dagger} + \cosh{\left (\chi \right )} {a}\right)^{2} + \left(\sinh{\left (\chi \right )} {a} + \cosh{\left (\chi \right )} {{a}^\dagger}\right)^{2}\right)$$
In [38]:
H6 = normal_ordered_form(H5.expand(), independent=True)

H6 = drop_c_number_terms(H6)

H6 = collect(H6, [Dagger(a) * a, Dagger(a)**2, a**2])

H6

Out[38]:
$$\left(\Delta_{a} \sinh^{2}{\left (\chi \right )} + \Delta_{a} \cosh^{2}{\left (\chi \right )} + 4 \kappa \sinh{\left (\chi \right )} \cosh{\left (\chi \right )}\right) {{a}^\dagger} {a} + \left(\Delta_{a} \sinh{\left (\chi \right )} \cosh{\left (\chi \right )} + \kappa \sinh^{2}{\left (\chi \right )} + \kappa \cosh^{2}{\left (\chi \right )}\right) \left({{a}^\dagger}\right)^{2} + \left(\Delta_{a} \sinh{\left (\chi \right )} \cosh{\left (\chi \right )} + \kappa \sinh^{2}{\left (\chi \right )} + \kappa \cosh^{2}{\left (\chi \right )}\right) \left({a}\right)^{2}$$
In [39]:
H7 = collect(H6, H6.args[1].args[0])

H7

Out[39]:
$$\left(\Delta_{a} \sinh^{2}{\left (\chi \right )} + \Delta_{a} \cosh^{2}{\left (\chi \right )} + 4 \kappa \sinh{\left (\chi \right )} \cosh{\left (\chi \right )}\right) {{a}^\dagger} {a} + \left(\Delta_{a} \sinh{\left (\chi \right )} \cosh{\left (\chi \right )} + \kappa \sinh^{2}{\left (\chi \right )} + \kappa \cosh^{2}{\left (\chi \right )}\right) \left(\left({{a}^\dagger}\right)^{2} + \left({a}\right)^{2}\right)$$
In [40]:
# Trick to simplify the coefficients for the quantum operators
H8 = Add(*(simplify(arg.args[0]) * Mul(*(arg.args[1:])) for arg in H7.args))

H8

Out[40]:
$$\left(\frac{\Delta_{a}}{2} \sinh{\left (2 \chi \right )} + \kappa \cosh{\left (2 \chi \right )}\right) \left(\left({{a}^\dagger}\right)^{2} + \left({a}\right)^{2}\right) + \left(\Delta_{a} \cosh{\left (2 \chi \right )} + 2 \kappa \sinh{\left (2 \chi \right )}\right) {{a}^\dagger} {a}$$

Now if we choose $\chi$ such that the coefficient of $a^2 + {a^\dagger}^2$ is zero:

In [41]:
chi_eq = H8.args[0].args[0]

In [42]:
Eq(chi_eq, 0)

Out[42]:
$$\frac{\Delta_{a}}{2} \sinh{\left (2 \chi \right )} + \kappa \cosh{\left (2 \chi \right )} = 0$$
In [43]:
chi_sol = simplify(solve(chi_eq, chi)[3])

Eq(chi, chi_sol)

Out[43]:
$$\chi = \log{\left (\sqrt[4]{\frac{\Delta_{a} - 2 \kappa}{\Delta_{a} + 2 \kappa}} \right )}$$

we obtain a diagonal Hamiltonian:

In [44]:
H9 = H8.subs(chi_eq.args[0], -chi_eq.args[1])

H9

Out[44]:
$$\left(\Delta_{a} \cosh{\left (2 \chi \right )} + 2 \kappa \sinh{\left (2 \chi \right )}\right) {{a}^\dagger} {a}$$

with frequency:

In [45]:
Eq(Symbol("omega"), H9.args[0])

Out[45]:
$$\omega = \Delta_{a} \cosh{\left (2 \chi \right )} + 2 \kappa \sinh{\left (2 \chi \right )}$$

### Versions¶

In [46]:
%reload_ext version_information

%version_information sympy, sympsi

Out[46]:
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
Thu Oct 09 16:01:18 2014 JST