Reproducing some results from K.D. Petersson et al. (Nature 2012)

In [1]:
%pylab inline
Populating the interactive namespace from numpy and matplotlib
In [2]:
from sympy import *
import numpy
In [3]:
init_printing()
In [4]:
TP = 2 * numpy.pi
In [5]:
ai = Symbol('a_in')
ao = Symbol('a_out')
bi = Symbol('b_in')
bo = Symbol('b_out')
In [6]:
# fundamental parameters
wR = Symbol('omega_R', real=True) # probe field frequency
w0 = Symbol('omega_0', real=True) # cavity frequency

k1 = Symbol('kappa_1', real=True)
k2 = Symbol('kappa_2', real=True)
ki = Symbol('kappa_i', real=True)
sigma_e = Symbol('sigma_e', real=True)

gamma  = Symbol('gamma', real=True)

tc = Symbol('t_c', real=True)
de = Symbol('delta_epsilon', real=True)
e  = Symbol('epsilon', real=True)
gc = Symbol('g_c', real=True)

X, Y = symbols("X, Y")
In [7]:
param_fundamental     = {tc: 3.5 * TP, 
                         w0: 6.1948 * TP,
                         gc: 0.030 * TP, 
                         k1: 0.015, 
                         k2: 0.015, 
                         ki: 0.00, 
                         sigma_e: 5.1 * TP, 
                         gamma: 0.066} # 15ns

param_fundamental_all = {wR: 6.1948 * TP,
                         e: 0}

param_fundamental_all.update(param_fundamental)
In [8]:
# derived parameters
delta0 = Symbol('Delta_0', real=True)
delta  = Symbol('Delta',   real=True)
W      = Symbol('Omega', real=True)

geff   = Symbol('g_eff', real=True)
chi    = Symbol('chi', real=True)
k      = Symbol('kappa', real=True)

Using derived parameters

In [9]:
ao = -I * sqrt(k1*k2) * bi / (delta0 - I*k/2 + geff * chi)

ao
Out[9]:
$$- \frac{i b_{in} \sqrt{\kappa_{1} \kappa_{2}}}{\Delta_{0} + \chi g_{eff} - \frac{i \kappa}{2}}$$
In [10]:
# the transmission coefficient
T = ao.subs(ai,0) / (bi)

T
Out[10]:
$$- \frac{i \sqrt{\kappa_{1} \kappa_{2}}}{\Delta_{0} + \chi g_{eff} - \frac{i \kappa}{2}}$$

Using fundamental parameters

In [11]:
W_num = sqrt((2*tc)**2 + e**2)

W_num
Out[11]:
$$\sqrt{\epsilon^{2} + 4 t_{c}^{2}}$$
In [12]:
# probe field detuning from atom
delta_num = W - wR

delta_num
Out[12]:
$$\Omega - \omega_{R}$$
In [13]:
# cavity frequency detuning from probe frequency
delta0_num = w0 - wR

delta0_num
Out[13]:
$$\omega_{0} - \omega_{R}$$
In [14]:
geff_num = gc * 2 * tc / W

geff_num
Out[14]:
$$\frac{2 g_{c}}{\Omega} t_{c}$$
In [15]:
chi_num = geff / (-delta + I * gamma / 2)

chi_num
Out[15]:
$$\frac{g_{eff}}{- \Delta + \frac{i \gamma}{2}}$$
In [16]:
delta_phi = - arg(T) * 180/pi
In [17]:
param_derived = {W: W_num,
                 geff: geff_num,
                 chi:  chi_num,
                 delta0: delta0_num,
                 delta: delta_num,
                 k: k1 + k2 + ki}
In [18]:
chi.subs(param_derived).subs(param_derived)
Out[18]:
$$\frac{2 g_{c} t_{c}}{\sqrt{\epsilon^{2} + 4 t_{c}^{2}} \left(\frac{i \gamma}{2} + \omega_{R} - \sqrt{\epsilon^{2} + 4 t_{c}^{2}}\right)}$$
In [19]:
T.subs(param_derived).subs(param_derived)
Out[19]:
$$- \frac{i \sqrt{\kappa_{1} \kappa_{2}}}{\frac{4 g_{c}^{2} t_{c}^{2}}{\left(\epsilon^{2} + 4 t_{c}^{2}\right) \left(\frac{i \gamma}{2} + \omega_{R} - \sqrt{\epsilon^{2} + 4 t_{c}^{2}}\right)} + \omega_{0} - \omega_{R} - \frac{i}{2} \left(\kappa_{1} + \kappa_{2} + \kappa_{i}\right)}$$
In [20]:
delta_phi.subs(param_derived).subs(param_derived)
Out[20]:
$$- \frac{180}{\pi} \arg{\left (- \frac{i \sqrt{\kappa_{1} \kappa_{2}}}{\frac{4 g_{c}^{2} t_{c}^{2}}{\left(\epsilon^{2} + 4 t_{c}^{2}\right) \left(\frac{i \gamma}{2} + \omega_{R} - \sqrt{\epsilon^{2} + 4 t_{c}^{2}}\right)} + \omega_{0} - \omega_{R} - \frac{i}{2} \left(\kappa_{1} + \kappa_{2} + \kappa_{i}\right)} \right )}$$

Test conditions

In [21]:
(2*tc > wR).subs(param_derived).subs(param_fundamental_all)
Out[21]:
True
In [22]:
(delta > gamma/2).subs(param_derived).subs(param_fundamental_all)
Out[22]:
True

Using numerical values for parameters

In [23]:
(2*tc/wR).subs(param_derived).subs(param_fundamental_all).evalf()
Out[23]:
$$1.12997998321173$$
In [24]:
(delta/(gamma/2)).subs(param_derived).subs(param_fundamental_all).evalf()
Out[24]:
$$153.309721495182$$

Unit conversions

In [25]:
def GHz_to_meV(w_GHz):
    # 1 GHz = 4.1357e-6 eV = 4.1357e-3 meV
    w_meV = w_GHz * 4.1357e-3
    return w_meV

def meV_to_GHz(w_meV):
    # 1 meV = 1.0/4.1357e-3 GHz
    w_GHz = w_meV / 4.1357e-3
    return w_GHz

Plot the reflection phase shift as a function of probe freqency

In [26]:
wR_vec = linspace(6.190, 6.20, 1000) * TP

T0_lambda = lambdify(X, T.subs(param_derived).subs(param_derived).subs({e: 0.0, wR: X}).subs(param_fundamental), 'numpy')
T1_lambda = lambdify(X, T.subs(param_derived).subs(param_derived).subs({e: 3.0 * TP, wR: X}).subs(param_fundamental), 'numpy')

T0_vec = T0_lambda(wR_vec)
T1_vec = T1_lambda(wR_vec)
In [27]:
fig, axes = subplots(1,2, figsize=(12,4))

# according to phase definition in Petersson et al
axes[0].plot(wR_vec/TP, -numpy.angle(T0_vec) * 180/pi, color="green")
axes[0].plot(wR_vec/TP, -numpy.angle(T1_vec) * 180/pi, 'b:')
axes[0].set_ylabel("phase", fontsize=18)
axes[0].set_xlabel(r'$\omega_R$', fontsize=18)

axes[1].plot(wR_vec/TP, numpy.abs(T0_vec)**2, color="green")
axes[1].plot(wR_vec/TP, numpy.abs(T1_vec)**2, 'b:')
axes[1].set_ylabel(r'$|T|^2$', fontsize=18)
axes[1].set_xlabel(r'$\omega_R$', fontsize=18);

fig.tight_layout()

Plot the reflection phase shift as a function of bias

In [28]:
e_vec = linspace(-20.0, 20.0, 500) * TP

T_lambda = lambdify(X, T.subs(param_derived).subs(param_derived).subs({tc: (3.5/2) * TP, e: X}).subs(param_fundamental_all), 'numpy')

T_vec = T_lambda(e_vec)
In [29]:
fig, axes = subplots(1,2, figsize=(12,4))

# according to phase definition in Petersson et al
axes[0].plot(e_vec/TP, -numpy.angle(T_vec) * 180/pi)
axes[1].plot(e_vec/TP, numpy.abs(T_vec)**2)

axes[0].set_ylabel("phase", fontsize=18)
axes[0].set_xlabel(r'$\epsilon$ (GHz)', fontsize=18)
axes[1].set_ylabel(r'$|T|^2$', fontsize=18)
axes[1].set_xlabel(r'$\epsilon$ (GHz)', fontsize=18);

fig.tight_layout()

Plot Omega vs bias

In [30]:
tc_vec = array([1.8, 5.1, 5.8, 6.2, 7.0]) / 2.0 * TP
In [31]:
e_vec = np.linspace(meV_to_GHz(-0.12), meV_to_GHz(0.12), 1000) * TP
In [32]:
W_lambda = lambdify((X, Y), W.subs(param_derived).subs(param_derived).subs(e, X).subs(tc, Y), 'numpy')
W_mat = [W_lambda(e_vec, tc_num) for tc_num in tc_vec]

tc_num = tc_vec[0]
chi_lambda = lambdify((X, Y), chi.subs(param_derived).subs(param_derived).subs({e: X, tc: Y}).subs(param_fundamental_all), 'numpy')
chi_vec = chi_lambda(e_vec, tc_num)
In [33]:
w0_num = (w0 / (2*pi)).subs(param_fundamental_all).evalf()
In [34]:
fig, axes = subplots(1,2, figsize=(10, 5))

axes[0].plot(GHz_to_meV(e_vec/TP), +0.5 * real(W_mat[0])/TP)
axes[0].plot(GHz_to_meV(e_vec/TP), -0.5 * real(W_mat[0])/TP)
axes[0].set_ylabel(r'$E_\pm$', fontsize=18)
axes[0].set_xlabel(r'$\epsilon$ (meV)', fontsize=18)

axes[1].plot(GHz_to_meV(e_vec/TP), real(chi_vec))
axes[1].set_ylabel(r'$\chi$', fontsize=18)
axes[1].set_xlabel(r'$\epsilon$ (meV)', fontsize=18)

fig.tight_layout()
In [35]:
fig, axes = subplots(1,2, figsize=(10,5))

# GHz
for W_vec in W_mat:
    axes[0].plot(e_vec/TP, real(W_vec)/TP)

axes[0].plot(e_vec/TP, w0_num * np.ones(shape(e_vec)), 'k--')

axes[0].set_ylabel(r'$\Omega$', fontsize=18)
axes[0].set_xlabel(r'$\epsilon$ (GHz)', fontsize=18)
axes[0].set_ylim(0, 11)
axes[1].set_xlim(meV_to_GHz(-0.05), meV_to_GHz(0.05))

# meV
for W_vec in W_mat:
    axes[1].plot(GHz_to_meV(e_vec/TP), real(W_vec)/TP)

axes[1].plot(GHz_to_meV(e_vec/TP), w0_num * np.ones(shape(e_vec)), 'k--')

axes[1].set_ylabel(r'$\Omega$', fontsize=18)
axes[1].set_xlabel(r'$\epsilon$ (meV)', fontsize=18)
axes[1].set_ylim(0, 11)
axes[1].set_xlim(-0.05, 0.05)

fig.tight_layout()

Plot phase shift vs bias

In [36]:
#delta_phi_mat = array([[delta_phi.subs(param_derived).subs({e: e_num, tc: tc_num}).subs(param_fundamental_all).evalf()
#                        for e_num in e_vec]
#                       for tc_num in tc_vec], dtype=complex)

delta_phi_lambda = lambdify((X, Y), delta_phi.subs(param_derived).subs(param_derived).subs({e: X, tc: Y}).subs(param_fundamental_all), 'numpy')
delta_phi_mat = array([delta_phi_lambda(e_vec, tc_num) for tc_num in tc_vec], dtype=complex)

Without cavity convolution

In [37]:
fig, axes = subplots(1,2, figsize=(10,5))

# GHz
for idx, delta_phi_vec in enumerate(delta_phi_mat):
    axes[0].plot(e_vec/TP, - idx * 0.0 + real(delta_phi_vec))

axes[0].set_ylabel(r'$\Delta\varphi$', fontsize=18)
axes[0].set_xlabel(r'$\epsilon$ (GHz)', fontsize=18)
axes[1].set_xlim(meV_to_GHz(-0.05), meV_to_GHz(0.05))

# meV
for idx, delta_phi_vec in enumerate(delta_phi_mat):
    axes[1].plot(GHz_to_meV(e_vec/TP), real(delta_phi_vec))

axes[1].set_ylabel(r'$\Delta\varphi$', fontsize=18)
axes[1].set_xlabel(r'$\epsilon$ (meV)', fontsize=18)
axes[1].set_xlim(-0.05, 0.05)

fig.tight_layout()

Definition of cavity convolution function

In [38]:
def cavity_convolution_func(w_vec, w0):
    g = numpy.exp(-(w_vec-w0)**2/param_fundamental[sigma_e]**2)

    return g / numpy.sum(g)
In [39]:
def cavity_convolution(w_vec, f_vec):
    return array([numpy.sum(f_vec * cavity_convolution_func(w_vec, w)) for w in w_vec])
In [40]:
fig, ax = subplots(1, 1, figsize=(10,5))

f_vec = abs(e_vec) < 5*TP
w0 = 0 * TP
ax.plot(e_vec/TP, cavity_convolution_func(e_vec, w0))
ax.plot(e_vec/TP, f_vec)
ax.plot(e_vec/TP, cavity_convolution(e_vec, f_vec));

With cavity convolution

In [41]:
fig, axes = subplots(1,2, figsize=(10,5))

# GHz
for idx, delta_phi_vec in enumerate(delta_phi_mat):
    axes[0].plot(e_vec/TP, - idx * 7.5 + real(cavity_convolution(e_vec, delta_phi_vec)))

axes[0].set_ylabel(r'$\Delta\varphi$', fontsize=18)
axes[0].set_xlabel(r'$\epsilon$ (GHz)', fontsize=18)
axes[1].set_xlim(meV_to_GHz(-0.05), meV_to_GHz(0.05))

# meV
for idx, delta_phi_vec in enumerate(delta_phi_mat):
    axes[1].plot(GHz_to_meV(e_vec/TP), - idx * 7.5 + real(cavity_convolution(e_vec, delta_phi_vec)))

axes[1].set_ylabel(r'$\Delta\varphi$', fontsize=18)
axes[1].set_xlabel(r'$\epsilon$ (meV)', fontsize=18)
axes[1].set_xlim(-0.05, 0.05)

fig.tight_layout()

Versions

In [42]:
%load_ext version_information
In [43]:
%version_information numpy, sympy, matplotlib
Out[43]:
SoftwareVersion
Python3.3.2+ (default, Oct 9 2013, 14:50:09) [GCC 4.8.1]
IPython1.1.0
OSposix [linux]
numpy1.7.1
sympy0.7.3
matplotlib1.3.1
Fri Nov 01 11:53:00 2013 KST