Yanbing Liu and Andrew A Houck, Quantum electrodynamics near a photonic bandgap, Nat. Phys., (2016).
Reproduced by Eunjong Kim (eunjongkim@caltech.edu)
# Setup for symbolic calculation
import sympy as sp
sp.init_printing()
# Setup for numerical calculation
import numpy as np
import numpy.linalg as la
from scipy.optimize import fsolve, brentq
from scipy.integrate import quad
# Setup for plots
import matplotlib.pyplot as plt
from matplotlib import rcParams, cm
%matplotlib inline
rcParams.update({"text.usetex": True, "font.size": 16})
In this section, we derive the equation (1) of this paper using the SymPy module (symbolic calulation).
# Define symbolic variables corresponding to frequency,
# wavevector, and size of the unit cell
omega_sym = sp.Symbol(r'\omega')
k_sym = sp.Symbol(r'k')
d_sym = sp.Symbol(r'd')
# Define a list of phase velocity, characteristic impedance,
# and length of each section of the unit cell
v_sym = [sp.Symbol(r'v_{p1}'), sp.Symbol(r'v_{p2}')]
Z_sym = [sp.Symbol(r'Z_{1}'), sp.Symbol(r'Z_{2}')]
l_sym = [sp.Symbol(r'l_{1}'), sp.Symbol(r'l_{2}')]
Z0_sym = sp.Symbol(r'Z_{0}')
# qubit parameters
gamma_sym = sp.Symbol(r'\gamma')
omega_a_sym = sp.Symbol(r'\omega_a')
We define the transfer matrix (ABCD) of each coplanar waveguide (CPW) section,
def CPW_ABCD_sym(omega, idx):
"""
Construct an ABCD matrix (symbolic) of the coplanar waveguide of specified index
"""
A = sp.cos(omega * l_sym[idx] / v_sym[idx])
B = sp.I * sp.sin(omega * l_sym[idx] / v_sym[idx]) * Z_sym[idx]
C = sp.I * sp.sin(omega * l_sym[idx] / v_sym[idx]) / Z_sym[idx]
D = sp.cos(omega * l_sym[idx] / v_sym[idx])
return sp.Matrix([[A, B], [C, D]])
The transfer matrix $M_\textrm{unit}$ is the product of that of each section:
M_unit_sym = CPW_ABCD_sym(omega_sym, 0) * CPW_ABCD_sym(omega_sym, 1)
M_unit_sym
Also, the ABCD matrix $M_a$ for qubit can be described as:
def qubit_ABCD_sym(omega):
"""
Construct an ABCD matrix (symbolic) of an atom in the single photon regime
"""
A = 1
B = 0
C = - sp.I * gamma_sym / (omega - omega_a_sym) / Z0_sym
D = 1
return sp.Matrix([[A, B], [C, D]])
qubit_ABCD_sym(omega_sym)
Now, we calculate the band structure for the bare photonic crystal (in the absence of qubit). For a periodic system, the voltage $V(x, t)$ and current $I(x, t)$ should be of the Bloch form, $$ \begin{bmatrix} V(0, t) \\ I(0, t) \end{bmatrix} = M_\textrm{unit} \begin{bmatrix} V(d, t) \\ I(d, t) \end{bmatrix} = M_\textrm{unit}e^{ikd}\begin{bmatrix} V(0, t) \\ I(0, t) \end{bmatrix}, $$ so that $$ M_\textrm{unit}\begin{bmatrix} V(0, t) \\ I(0, t) \end{bmatrix} = e^{-ikd}\begin{bmatrix} V(0, t) \\ I(0, t) \end{bmatrix}. $$
We require that the determinant of $ (e^{ikd} M_\textrm{unit} - I)$ be zero for non-trivial solutions, and we get the following dispersion relation:
Expr = sp.simplify(sp.det(M_unit_sym * sp.exp(sp.I * k_sym * d_sym) - sp.eye(2))
/ (2 * sp.exp(sp.I * k_sym * d_sym))).subs(d_sym, sum(l_sym))
sp.Eq(Expr, 0)
Assuming the same phase velocity in all sections, $v_{p} \equiv v_{p1} = v_{p2}$
v_p_sym = sp.Symbol('v_p')
Expr = Expr.subs({v_sym[0]: v_p_sym, v_sym[1]: v_p_sym})
sp.Eq(Expr, 0)
Thus, the Equation (1) is obtained.
From this dispersion relation, it is possible to numerically evaluate the band structure. First, we put the values of parameters. The eigenfrequency of the second band at $k=\pi/d$ is given by $\omega_0 /2\pi= 7.7\textrm{ GHz}$. The parameters corresponding to lengths $l$ and the impedances $Z$ are: $l_1 = 0.45\mathrm{\ mm}$, $l_2 = 8\mathrm{\ mm}$, $Z_1 = 28\mathrm{\ \Omega}$, $Z_2 = 125\mathrm{\ \Omega}$.
# values of parameters in SI units
omega0 = 2 * np.pi * 7.7e9
l = [0.45e-3, 8e-3]
Z = [28, 125]
d = sum(l)
dict_items = []
for idx in range(2):
dict_items.append((l_sym[idx], l[idx]))
dict_items.append((Z_sym[idx], Z[idx]))
We can compute the value of the phase velocity $v_p$ (assuming $v_{p1} = v_{p2}$) using the dispersion relation and $\omega = \omega_0$ at $k = \pi/d$
Expr_ = Expr.subs(dict_items)
expr_ = Expr_.subs({omega_sym: omega0, k_sym: np.pi/d}); expr_
func = sp.lambdify(v_p_sym, expr_, "numpy")
v_p = fsolve(func, omega0/(np.pi/d))[0]; v_p
Thus the value of the phase velocity is estimated as $v_p = 1.248 \times 10^{8} \mathrm{\ m/s}$. Using this, we calculate the band structure from the dispersion relation
expr = Expr_.subs(v_p_sym, v_p)
func = sp.lambdify((omega_sym, k_sym), expr, "numpy")
# Number of k-points to evaluate
N_kPts = 100
# Number of bands to evaluate
N_band = 4
kPts = np.linspace(- np.pi / d, np.pi / d, N_kPts)
freq = np.zeros((N_band, N_kPts))
for n in range(N_band):
for idx, k_ in enumerate(kPts):
# search for eigenfrequency between
# n v_p k_0 and (n+1) v_0 k_0 (k_0 = \pi / d)
freq[n, idx] = brentq(func, n * (v_p * np.pi / d),
(n + 1) * (v_p * np.pi / d),
args=(k_)) / (2 * np.pi)
fig, ax = plt.subplots(1, 1, figsize=(5, 4))
fig_colors = ['green','red','blue', 'purple']
for n in range(N_band):
ax.plot(kPts / (np.pi / d), freq[n, :] / 1e9, color=fig_colors[n])
ax.set_xlabel(r'Wavevector $k/(\pi/d)$');
ax.set_ylabel('Frequency [GHz]');
# characteristic impedance of the measurement chain
Z0 = 50
def ABCDtoS(ABCD):
A, B = ABCD[0]
C, D = ABCD[1]
denom = (A + B / Z0 + C * Z0 + D)
S11 = (A + B / Z0 - C * Z0 - D) / denom
S12 = 2 * (A * D - B * C) / denom
S21 = 2 / denom
S22 = (- A + B / Z0 - C * Z0 + D) / denom
return np.array([[S11, S12], [S21, S22]])
Here, I calculate the transmission amplitude of the bare photonic crystal (14 cells) using the transfer matrix method.
M_unit_sym_ = M_unit_sym.subs(dict_items).subs({v_sym[0]: v_p, v_sym[1]: v_p})
M_unit = sp.lambdify(omega_sym, M_unit_sym_, "numpy")
N_cell = 14 # Number of unit cells to calculate
N_freqPts = 100
freqPts = np.linspace(5.5, 8.5, N_freqPts) # frequency points in units of GHz
S12dB = np.zeros(N_freqPts)
for idx, freq_ in enumerate(freqPts):
S12dB[idx] = 20 * np.log10(ABCDtoS(la.matrix_power(
M_unit(2e9*np.pi*freq_), N_cell))[0, 1])
/opt/local/Library/Frameworks/Python.framework/Versions/3.4/lib/python3.4/site-packages/ipykernel/__main__.py:8: ComplexWarning: Casting complex values to real discards the imaginary part
fig, ax = plt.subplots(1, 1, figsize=(5, 3))
ax.plot(freqPts, S12dB, color='red', lw=2)
ax.set_xlabel('Frequency [GHz]')
ax.set_ylabel('$S_{12}$ [dB]')
ax.set_ylim([-40, 5]);
To analyze the mode structure (Bloch wavefunction) of an infinite microwave photonic crystal, we solve the following series of differential equations
$$ \frac{\partial}{\partial x} V(x, t) = -l(x) \frac{\partial}{\partial t} I(x,t), \quad \frac{\partial}{\partial x} I(x, t) = -c(x) \frac{\partial}{\partial t} V(x,t) $$which reduces to a single partial differential equation of $V(x, t)$
$$ \frac{\partial}{\partial x} \left[ \frac{v_p}{Z_c(x)} \frac{\partial V(x,t)}{\partial x}\right] = \frac{1}{v_p Z_c(x)} \frac{\partial^2}{\partial t^2} V(x,t) $$The characteristic impedance $Z_c(x) = \sqrt{l(x)/c(x)}$ inside the unit cell ($-d/2< x < d/2$) is let to be $$ Z_c (x) = $$ while the phase velocity $v_p = 1/\sqrt{l(x)c(x)}$ is a constant.
The time-independent ODE is obtained by separation of variables $V(x, t) = V(x) e^{i\omega t}$ $$ \frac{d}{dx} \left[\frac{1}{Z_c(x)} \frac{dV(x)}{dx}\right] = - \left(\frac{\omega}{v_p}\right)^2 \frac{V(x)}{Z_c(x)} $$
We solve the partial differential equation
For numerical calculation of the mode structure, we impose the periodic boundary condition $Z_c(x) = Z_c(x+d)$ Fourier series expand $1/Z_c(x)$ as $$ \frac{1}{Z_c (x)} = \sum_{m} \eta_{m} e^{i2\pi m x/d}. $$
Here the Fourier coefficients $\eta_m$ are written as $$\eta_m = \frac{1}{d} \int_{-d/2}^{d/2} dx \frac{e^{-i 2\pi m x/d}}{Z_c(x)}$$
# Characteristic impedance inside the unit cell
def Zc(x):
return Z[1] * ((x > - l[1] / 2) * (x < l[1] / 2)) \
+ Z[0] * ((x < - l[1] / 2) + (x > l[1] / 2))
def eta(n_):
"""
The n-th discrete Fourier coefficient of 1/Zc(x).
Imaginary part vanishes due to the even parity of Zc(x)
with respect to the x = 0 plane.
"""
integrand = lambda x: (1 / Zc(x) *
np.exp(- 2j * np.pi * n_ * x / d))
return quad(integrand, - d / 2, d / 2)[0] / d
N_xPts = 100
xPts = np.linspace(- d / 2, d / 2, N_xPts)
# discrete Fourier space dimension
L = 10
n = np.arange(- 2 * L, 2 * L + 1)
#
eta_array = np.zeros(4 * L + 1)
for n_ in n:
eta_array[n_ + 2 * L] = eta(n_)
/opt/local/Library/Frameworks/Python.framework/Versions/3.4/lib/python3.4/site-packages/scipy/integrate/quadpack.py:380: ComplexWarning: Casting complex values to real discards the imaginary part return _quadpack._qagse(func,a,b,args,full_output,epsabs,epsrel,limit)
fig, axes = plt.subplots(1, 2, figsize=(10, 4))
axes[0].plot(xPts / d, Zc(xPts), color=fig_colors[0], lw=2)
axes[0].set_xlabel('Position in unit cell $x/d$')
axes[0].set_ylabel(r'Characteristic Impedance $Z_c$ [$\Omega$]')
axes[0].set_xlim([- 1/2, 1/2]);
axes[1].bar(n, eta_array, color=fig_colors[2]);
axes[1].set_xlabel(r'$m$')
axes[1].set_xlim([- 2 * L, 2 * L]);
axes[1].set_ylabel(r'Discrete Fourier Coefficient $\eta_m$ [$\Omega^{-1}$]')
fig.tight_layout();
To solve this differential equation, we introduce eigenfunctions of the Bloch form $V_k(x) = u_k(x) e^{ikx}$ which satisfies $u_k(x+d) = u_k(x)$ and thus can be Fourier series expanded as $$ u_k(x) = \sum_{m} c_{k,m} e^{i2\pi m x/d}$$
We introduce a new index $n$ to replace $m = n-m'$ $$ -\sum_{m', n} \eta_{n-m'} \left(k+\frac{2\pi m'}{d}\right)\left(k+\frac{2\pi n}{d}\right) c_{k,m'} e^{i(k+2\pi n /d)x} = - \left(\frac{\omega}{v_p}\right)^2 \sum_{m',n} \eta_{n-m'} c_{k,m'} e^{i(k+ 2\pi n /d)x} $$
Applying the inverse Fourier transform $\frac{1}{d}\int_{-d/2}^{d/2} dx \ e^{-i(k+2\pi m/d)x}$, $$ -\sum_{m', n} \eta_{n-m'} \left(k+\frac{2\pi m'}{d}\right)\left(k+\frac{2\pi n}{d}\right) c_{k,m'}\left[ \frac{1}{d}\int_{-d/2}^{d/2}dx\ e^{i2\pi (n-m)x /d }\right] = - \left(\frac{\omega}{v_p}\right)^2 \sum_{m',n} \eta_{n-m'} c_{k,m'}\left[ \frac{1}{d}\int_{-d/2}^{d/2}dx\ e^{i2\pi (n-m)x /d }\right] $$
The delta function identity $$ \frac{1}{d}\int_{-d/2}^{d/2} dx\ e^{i2\pi (n-m)x /d} = \delta_{n,m} $$ reduces the equation into a summation only in one index $m'$ $$ \sum_{m'} \eta_{m-m'} \left(k + \frac{2\pi m}{d}\right)\left(k + \frac{2\pi m'}{d}\right) c_{k, m'} = \left(\frac{\omega}{v_p}\right)^2 \sum_{m'} \eta_{m-m'} c_{k,m'} $$
By identifying the matrices $\mathbf{A}$, $\mathbf{B}$, and $\mathbf{c}$ as
$$[\mathbf{A}]_{m, m'} = \eta_{m-m'},\quad [\mathbf{B}]_{m, m'} = \eta_{m-m'} \left(k + \frac{2\pi m}{d}\right) \left(k + \frac{2\pi m'}{d}\right),\quad [\mathbf{c}]_m' = c_{k, m'},$$This is a matrix equation $\mathbf{B}\mathbf{c} = (\omega/v_p)^2 \mathbf{A}\mathbf{c}$ which can be viewed as an eigenequation $$(\mathbf{A}^{-1}\mathbf{B})\mathbf{c} = \left(\frac{\omega}{v_p}\right)^2 \mathbf{c}$$
A = np.zeros((2 * L + 1, 2 * L + 1))
for i in range(2 * L + 1):
for j in range(2 * L + 1):
m, m_ = i - L, j - L # index offset
A[i, j] = eta_array[m - m_ + 2 * L]
def B(k):
temp = np.zeros((2 * L + 1, 2 * L + 1))
for i in range(2 * L + 1):
for j in range(2 * L + 1):
m, m_ = i - L, j - L # index offset
temp[i, j] = (eta_array[m - m_ + 2 * L]
* (k + 2 * np.pi * m / d)
* (k + 2 * np.pi * m_ / d))
return temp
def Op(k):
return np.dot(la.inv(A), B(k))
def solve_eigenequation(k):
eig_vals, eig_vecs = la.eig(Op(k))
# sorting eigenvalues and eigenvectors
idx = eig_vals.argsort()
eig_vals = eig_vals[idx]
eig_vecs = eig_vecs[:,idx]
eig_freq = np.sqrt(eig_vals) * v_p / (2 * np.pi)
return eig_freq, eig_vecs
def Bloch_wavefunc(x, c):
"""`1
Given the amplitudes [c],
Calculate the Bloch wavefunction
u_k (x) = sum_{n} c_{k, n} exp(i2\pi n x/d)
and return its absolute value.
"""
V = np.zeros(len(x))
for n in range(2 * L + 1):
V = V + c[n] * np.exp(2j * np.pi * (n - L) * x / d)
return abs(V)
X, Y = np.meshgrid(kPts, xPts)
N_modes = 2 # Number of modes to consider
eig_funcs = []
eig_freqs = []
for n in range(N_modes):
eig_funcs.append(np.zeros((len(kPts), len(xPts))))
eig_freqs.append(np.zeros(len(kPts)))
for idx, k_ in enumerate(kPts):
eig_vals, eig_vecs = solve_eigenequation(k_)
for n in range(N_modes):
eig_funcs[n][:, idx] = Bloch_wavefunc(xPts, eig_vecs[:, n])
eig_freqs[n][idx] = eig_vals[n]
fig, axes = plt.subplots(2, 2, figsize=(11, 8))
# Band structure
for n in range(N_modes):
axes[0, 0].plot(kPts/(np.pi/d), eig_freqs[n]/1e9, label=r'$n=%d$' % n, lw=2)
axes[0, 0].set_xlabel(r'Wavevector $k/(\pi/d)$')
axes[0, 0].set_ylabel(r'Frequency $\omega/2\pi$ [GHz]')
axes[0, 0].set_title('Band Structure')
axes[0, 0].legend(loc=0)
# Transmission Amplitude
axes[0, 1].plot(freqPts, S12dB, color='red', lw=2)
axes[0, 1].set_xlabel('Frequency [GHz]')
axes[0, 1].set_ylabel('$S_{12}$ [dB]')
axes[0, 1].set_ylim([-40, 5])
axes[0, 1].set_title('Transmission Amplitude')
# Mode structure
p0 = axes[1, 0].pcolor(X/(np.pi/d), Y/d, eig_funcs[0], cmap=cm.jet, vmin=0, vmax=1.5)
p1 = axes[1, 1].pcolor(X/(np.pi/d), Y/d, eig_funcs[1], cmap=cm.jet, vmin=0, vmax=1.5)
axes[1, 0].set_xlabel(r'Wavevector $k/(\pi/d)$'); axes[1, 1].set_xlabel(r'Wavevector $k/(\pi/d)$')
axes[1, 0].set_ylabel(r'Position in unit cell $x/d$')
axes[1, 0].set_xlim([-1, 1]); axes[1, 1].set_xlim([-1, 1])
axes[1, 0].set_ylim([-0.5, 0.5]); axes[1, 1].set_ylim([-0.5, 0.5])
axes[1, 0].set_title(r'Mode structure ($n=0$ band)')
axes[1, 1].set_title(r'Mode structure ($n=1$ band)');
fig.tight_layout()
fig.subplots_adjust(right=0.85)
cbar_ax = fig.add_axes([0.9, 0.15, 0.02, 0.3])
fig.colorbar(p1, cax=cbar_ax);
All panels of Figure S1 has been reproduced.
M_a_sym = qubit_ABCD_sym(omega_sym).subs({gamma_sym: 2e9 * np.pi * 0.5, Z0_sym: 28})
M_a = sp.lambdify((omega_a_sym, omega_sym), M_a_sym, "numpy")
omega_a = 2e9 * np.pi * 8.5
freqPts = np.linspace(7.2, 8.9, N_freqPts) # Frequency Points in units of GHz
S12dB = np.zeros(N_freqPts)
for idx, freq in enumerate(freqPts):
omega_ = 2e9 * np.pi * freq
M_PhC = la.matrix_power(M_unit(omega_), 7)
M = np.dot(M_PhC, np.dot(M_a(omega_a, omega_), M_PhC))
S12dB[idx] = 20 * np.log10(ABCDtoS(M)[0][1])
/opt/local/Library/Frameworks/Python.framework/Versions/3.4/lib/python3.4/site-packages/ipykernel/__main__.py:9: ComplexWarning: Casting complex values to real discards the imaginary part
fig, ax = plt.subplots(1, 1, figsize=(3, 5))
ax.plot(freqPts, S12dB, lw=2, color='blue')
ax.set_xlabel(r'Frequency [GHz]')
ax.set_ylabel(r'$S_{12}$ [dB]');
%reload_ext version_information
%version_information sympy, numpy, scipy, matplotlib
Software | Version |
---|---|
Python | 3.4.5 64bit [GCC 4.2.1 Compatible Apple LLVM 7.0.2 (clang-700.1.81)] |
IPython | 5.0.0 |
OS | Darwin 15.6.0 x86_64 i386 64bit |
sympy | 1.0 |
numpy | 1.11.1 |
scipy | 0.18.0 |
matplotlib | 1.5.2 |
Tue Aug 23 00:33:34 2016 PDT |