Reproduced by: Eunjong Kim (ekim7206@gmail.com)
In this notebook, I find the expression for the normal mode splitting of optomechanical system, following the supplementary information of [Gröblacher et. al., Nature 460, 724 (2009).]
from sympy import *
init_printing(use_latex='png')
from sympsi.operator import Operator
from sympsi.boson import BosonOp
from sympsi.dagger import Dagger
# Constants
d = symbols("Delta", real=True) # detuning
hbar, g, O, wp, wm = symbols("hbar, g, Omega, omega_+, omega_-", positive=True)
The linearized Hamiltonian for a driven cavity mode coupled via radiation pressure to a harmonically bound mirror is
$$ H = -\frac{\hbar\Delta}{2} (X_c^2 + P_c^2 ) + \frac{\hbar\Omega}{2} (X_m^2 + P_m^2 ) - \hbar g X_c X_m $$where $\Omega$ is the frequency of the mechanical resonator, $\Delta\equiv \omega_L - (\omega_c-\delta_{rp})$ is the effective detuning between the driving field and the cavity frequency$.($\delta_{rp}$ is the mean shift of the cavity frequency arising from radiation frequency.)
This can either be expressed in terms of the matrix equation with vector $\vec{R}\equiv [X_c,\ P_c,\ X_m,\ P_m]^T$:
$$ H =\frac{\hbar}{2} \vec{R}^T M \vec{R} $$where $M$ is the matrix defined as
M = Matrix([[-d, 0, -g, 0],
[0, -d, 0, 0],
[-g, 0, O, 0],
[0, 0, 0, O]])
M
# Operators
Xc, Xm = Operator('X_c'), Operator('X_m') # Position Quadratures
Pc, Pm = Operator('P_c'), Operator('P_m') # Momentum Quadratures
ac, am = BosonOp('a_c'), BosonOp('a_m') # Boson Operators
The quadrature operators are related to the annihilation and creation operators by the following linear equation: $$ \vec{R} = Q \vec{a}$$ where $\vec{a} \equiv [a_c,\ a^\dagger_c,\ a_m,\ a^\dagger_m]^T$, and $Q$ is the coefficient matrix
R = Matrix([Xc, Pc, Xm, Pm])
a = Matrix([ac, Dagger(ac), am, Dagger(am)])
Q = Matrix([[1, 1, 0, 0],
[-I, I, 0, 0],
[0, 0, 1, 1],
[0, 0, -I, I]])/sqrt(2)
Eq(R, MatMul(Q, a))
i.e., $$X_{\alpha} = \frac{a_{\alpha}+a^\dagger_{\alpha}}{\sqrt{2}},\quad P_{\alpha} = -i\frac{a_{\alpha}-a^\dagger_{\alpha}}{\sqrt{2}}\quad (\alpha = c, m)$$
The quadratures satisfy the canonical commutation relations $$[X_\alpha, P_\beta] = i \delta_{\alpha, \beta},\quad [X_\alpha, X_\beta] = [P_\alpha, P_\beta] = 0.$$ which are condensed into one simple equation $$ [\vec{R}_{i}, \vec{R}_{j}] = i J_{ij} $$ where $J$ is the matrix defined as
# Constraint representative of the canonical commutation relation
J = Matrix([[0, 1, 0, 0],
[-1, 0, 0, 0],
[0, 0, 0, 1],
[0, 0, -1, 0]])
J
Note that the matrix $J$ has simple properties: $\quad J^T = J^{-1} = -J, \quad J^2=-I$.
Our goal is to find a linear transformation $S$ that converts the original coordinates $\vec{R}$ into the normal mode coordinates $\vec{R}^{NM} = [X_+, P_+, X_-, P_-]^T = S\vec{R}$, whose Hamiltonian representation is $$ H = \frac{\hbar\omega_+}{2} \left( X_+^2 + P_+^2 \right) + \frac{\hbar \omega_-}{2} \left( X_-^2 + P_-^2\right) $$
# Operators in the normal mode coordinates
X_p, X_m = Operator('X_+'), Operator('X_-') # Position Quadratures
P_p, P_m = Operator('P_+'), Operator('P_-') # Momentum Quadratures
a_p, a_m = BosonOp('a_+'), BosonOp('a_-') # Boson Operators
Rnm = Matrix([X_p, P_p, X_m, P_m])
anm = Matrix([a_p, Dagger(a_p), a_m, Dagger(a_m)])
Eq(Rnm, MatMul(Q, anm))
This can either be expressed in terms of the matrix equation with vector $\vec{R}^{NM}$:
$$ H = \frac{\hbar}{2} \left(\vec{R}^{NM}\right)^T D \vec{R}^{NM} = \frac{\hbar}{2} \left(S\vec{R}\right)^T D \left(S\vec{R}\right) = \frac{\hbar}{2} \vec{R}^T \left(S^T D S\right)\vec{R} $$where $D$ is the diagonal matrix defined as
D = diag(wp, wp, wm, wm)
D
(NOTE: This transformation betweeen $M$ and $D$ is not in general a similarity transformation, so that $M$ and $D$ doesn't need not share a same set of eigenvalues.)
The latter property, which guarantees that the canonical commutation relations are conserved under the transformation, are often termed as a symplectic transformation.
The conditions 1) and 2) described above themselves doesn't say anything about the normal mode frequencies $\omega_+,\ \omega_-$ or the transformation matrix $S$. But, combining two matrix equations, we can get a meaningful relation useful in getting normal mode frequencies.
First, we can readily notice that the matrix $S$ should be invertible. This is due to the fact that $\det{J} = 1$, so the determinant of the rhs of condition 2),
$$\det{(SJS^T)} = \det{S} \cdot\det{J}\cdot \det{S^T} = \det{J} (\det{S})^2$$should also be 1. If $\det{S} = 0$, this cannot be the case. It naturally follows that $|\det{S}|=1$.
Second, multiplying $J^{-1} S^{-1}$ on the left of the both sides of 2), we get
$$ J^{-1} S ^{-1}J = -J S ^{-1} J = (-J S^{-1}) SJS^T = -J^2 S^T = S^T. $$Finally, substituting into 1),
$$ M = S^T D S = (J^{-1} S^{-1} J)DS = J^{-1} S^{-1}(JD) S \quad \Rightarrow \quad JM = S^{-1}(JD) S. $$we find that $JM$ and $JD$ are related by the similarity transformation. As a consequence, the matrices $JD$ and $JM$ share the same set of eigenvalues.
Now, let $$ A^{-1} (JM) A = B^{-1} (JD) B = \mathrm{diag}{(\lambda_1, \lambda_2, \lambda_3, \lambda_4)}. $$
(with this definition, $S = BA^{-1}$).
Extracting the eigenvectors and eigenvalues for two matrices $JM$ and $JD$,
# extracting eigenvectors and eigenvalues for J*M
[(e1,_,[v1]), (e2,_,[v2]), (e3,_,[v3]), (e4,_,[v4])] = (J * M).eigenvects()
e1, e2, e3, e4
# extracting eigenvectors and eigenvalues for J*D
[(f1,_,[w1]), (f2,_,[w2]), (f3,_,[w3]), (f4,_,[w4])] = (J * D).eigenvects()
f1, f2, f3, f4
there is a one-to-one correspondence between eigenvalues, which gives us the normal mode frequencies $\omega_{\pm}$:
wpsol = solve(e2-f2, wp)[0]
Eq(wp, wpsol)
wmsol = solve(e4-f4, wm)[0]
Eq(wm, wmsol)
# simpler expression without an imaginary number
wpsol = sqrt((d**2 + O**2 + sqrt((d**2 -O**2)**2 -4*d*O*g**2))/2)
wmsol = sqrt((d**2 + O**2 - sqrt((d**2 -O**2)**2 -4*d*O*g**2))/2)
Eq(wp, wpsol), Eq(wm, wmsol)
Now, try to find the transformation matrix $S$. First construct matrices $A$ and $B$ from the eigenvectors multiplied by arbitrary constants (which will be determined later)
# arbitrary constants
aa, bb, cc, dd = symbols("alpha, beta, gamma, delta")
x, y, z, w = symbols("x, y, z, w")
ee2, ff2 = simplify(sqrt(2)*e2), sqrt(2)*f2
ee4, ff4 = simplify(sqrt(2)*e4), sqrt(2)*f4
subs_dic = {ee2:ff2, ee2**2+2*d**2:ff2**2+2*d**2, ee4:ff4, ee4**2+2*d**2: ff4**2+2*d**2}
u1, u2, u3, u4 = (v1.subs(subs_dic), v2.subs(subs_dic),
v3.subs(subs_dic), v4.subs(subs_dic))
A = (x*u1).row_join(y*u2).row_join(z*u3).row_join(w*u4)
A = simplify(A)
A
B = (aa*w1).row_join(bb*w2).row_join(cc*w3).row_join(dd*w4)
B
and use the relation $S = B A^{-1}$.
Ainv = A.inv()
Ainv = simplify(Ainv)
S = simplify(B * Ainv)
S = Matrix([[simplify(S[i, j].factor()) for j in range(4)] for i in range(4)])
S
The aritrary constants can be determined by the condition 2), $J = SJS^T$. Calculating $SJS^T$, we get
SJST = simplify(S * J * S.T)
SJST = Matrix([[simplify(SJST[i,j].subs({wp:wpsol, wm:wmsol}).factor()) for j in range(4)] for i in range(4)])
The non-vanishing terms are (0, 1), (1, 0), (2, 3), (3, 2)-th elements and should be set to be equal to
Eq(1, simplify(SJST[0, 1].subs(subs_dic).factor()))
Eq(1, simplify(SJST[2, 3].subs(subs_dic).factor()))
I will choose the arbitrary constants to be
\begin{align} x &= y = \sqrt{\frac{\omega_+ (\omega_+^2 - \Delta^2)}{\Omega (\omega_+^2 - \omega_-^2)}},\\ z &= w = \sqrt{\frac{\omega_- (\Delta^2 - \omega_-^2)}{\Omega (\omega_+^2 - \omega_-^2)}},\\ \alpha&=\beta=\gamma=\delta = 1 \end{align}then, the transformation matrix $S$ becomes
S_ = S.subs({x: sqrt(wp * (wp**2 - d**2)/(O*(wp**2-wm**2))),
y: sqrt(wp * (wp**2 - d**2)/(O*(wp**2-wm**2))),
z: sqrt(wm * (d**2 - wm**2)/(O*(wp**2-wm**2))),
w: sqrt(wm * (d**2 - wm**2)/(O*(wp**2-wm**2))),
aa:1, bb:1, cc:1, dd:1})
S_ = Matrix([[simplify(S_[i,j].factor()).factor() for j in range(4)] for i in range(4)])
S_
Let's check if this matrix satisfies the two conditions we discussed above:
#check1
STDS_ = simplify(S_.T * D * S_)
STDS_ = Matrix([[simplify(STDS_[i,j].subs({wp:wpsol, wm:wmsol}).factor()) for j in range(4)] for i in range(4)])
STDS_
#check2
SJST_ = simplify(S_ * J * S_.T)
SJST_ = Matrix([[simplify(SJST_[i,j].subs({wp:wpsol, wm:wmsol}).factor()) for j in range(4)] for i in range(4)])
SJST_
Both conditions are fully satisfied with this transformation $S$.
The transformation between bosonic operators in different frames is readily obtained: $$ Q\vec{a}^{NM} = S \left(Q\vec{a}\right) \quad \Rightarrow \quad \vec{a}^{NM} = (Q^{-1} S Q) \vec{a} = T\vec{a} $$
T = simplify((Q.inv() * S_ * Q))
lhs = anm[0]
rhs = simplify((T * a)[0])
Eq(lhs, rhs)
lhs = anm[1]
rhs = simplify((T * a)[1])
Eq(lhs, rhs)
lhs = anm[2]
rhs = simplify((T * a)[2])
Eq(lhs, rhs)
lhs = anm[3]
rhs = simplify((T * a)[3])
Eq(lhs, rhs)
Here, I will consider the special case of the optomechanical system, where the driving field is tuned to the upper sideband ($\Delta = -\Omega$). The transformation matrix $S$ becomes
S__ = simplify(S_.subs({d: -O, wp: wpsol.subs(d, -O), wm: wmsol.subs(d, -O)}))
lhs = Rnm
rhs = MatMul(S__, R)
Eq(lhs, rhs)
so that
\begin{align} X_{\pm} &= \sqrt[4]{\frac{\Omega\pm g}{\Omega}} \frac{\mp X_c + X_p}{\sqrt{2}},\\ P_{\pm} &= \sqrt[4]{\frac{\Omega}{(\Omega\pm g)}} \frac{\mp P_c + P_p}{\sqrt{2}}. \end{align}It can be easily shown that this transformation is in consistent with the original Hamiltonian.
H = hbar * wp/2 *(X_p**2 + P_p**2) + hbar * wm/2 * (X_m**2 + P_m**2)
H
By substituting into the normal mode Hamiltonian,
H_ = H.subs({wp: wpsol.subs(d, -O), wm: wmsol.subs(d, -O),
X_p: (rhs.doit())[0],
P_p: (rhs.doit())[1],
X_m: (rhs.doit())[2],
P_m: (rhs.doit())[3],
})
simplify(H_.expand().factor()).subs({sqrt(((O+g)**2).expand()): O+g,
Xm * Xc: Xc * Xm,
Pm * Pc: Pc * Pm}).expand()
%reload_ext version_information
%version_information sympy, sympsi
Software | Version |
---|---|
Python | 3.4.1 (default, Sep 20 2014, 19:44:17) [GCC 4.2.1 Compatible Apple LLVM 5.1 (clang-503.0.40)] |
IPython | 2.3.0 |
OS | posix [darwin] |
sympy | 0.7.5-git |
sympsi | 0.1.0.dev-11eaf6c |
Fri Oct 24 15:05:38 2014 JST |