# Bose-Hubbard model in CQED systems

In this note, I will show how to realize effective Bose-Hubbard model in CQED systems. This note is mainly based on the paper Nature Physics 2, 849 – 855 (2006).

1. Introduction

Perturbative method can be well described the weak interation systems. But it is not useful in strong correlated systems. Besides, some key phenomenas in strong correlated systems are hardly ovserved because the time scale associated with them is very short. This is also the obstacles for realizing quantum processes.

Several important processes have been archieved in cold atomic systems, such as creating effective Bose-Hubbard model and realizing quantum phase transition in optical lattice trapped cold atoms. In this systems, a precise control of collective parameters is feasible. But it is still hard to indiviually address and control on single sites. Here a different approach to creat effective Bose-Hubbard model is proposed.

We consider an array of cavities. The dynamics of polaritons, combined atom photon excitations, in this setup is studied. Since the distance between adjucent cavities is considerably larger than the optial wavelength of the resonant mode., individual cavities can be addressed.

Model is shown bellow (refer to Fig 1 of quant-ph/0606097)

Photon hopping occurs between neighboring cavities. The repulsive force betweeen two polaritons occupying the same site is generated by a large Kerr nonlinearity. Three species of polaritons appear in our system. The species which is most steable is gonverned by effective Bose-Hubbard Hamiltonian.

$\displaystyle H_{eff} = \kappa \sum_{\vec{R}} (P_{\vec{R}}^\dagger)^2 (P_{\vec{R}})^2 + J \sum_{<\vec{R},\vec{R'}>} (P^\dagger_{\vec{R}} P_{\vec{R}'} + \mathrm{H.C.}), \ \ \ \ \ (1)$

where ${P^\dagger_{\vec{R}}}$ creates a polariton in the cavity of site ${\vec{R}}$, and the parameters ${\kappa}$ and ${J}$ describe on site repulsion and inter cavity hopping respectively.

2. Quantised field in a periodic array of cavities

In periodic array of cavities, we have periodic dielectric constant

$\displaystyle \epsilon(\vec{r}) = \epsilon(\vec{r}+ \vec{R}), \vec{R} = \vec{n} R, \ \ \ \ \ (2)$

where ${\epsilon}$ is real. Asorption is neglected. The field is represented by vector potential ${\vec{A}}$ and a scalar potential ${\Phi}$. They obey gauge condition: ${\Phi=0}$ and ${\Delta\cdot (\epsilon (\vec{r}) \vec{A}) = 0}$. Maxwell equation is

\displaystyle \begin{aligned} \vec{B}= \nabla \times \vec{A}; \quad \vec{E} = - \nabla \Phi - \frac{1}{c} \frac{\partial \vec{A}}{\partial t};\\ \vec{B}=\vec{H}; \quad \frac{1}{c} \frac{\partial \vec{D}}{\partial t} = \nabla \times \vec{H}. \end{aligned} \ \ \ \ \ (3)

Let’s set ${\vec{D}= \epsilon(\vec{r}) \vec{E}}$. Because of gauge choice, we get ${\vec{E}= -\frac{1}{c} \frac{\partial \vec{A}}{\partial t}}$. Then we get

$\displaystyle \frac{\epsilon(\vec{r})}{c^2} \frac{\partial^2 \vec{A}}{\partial t^2} + \nabla \times (\nabla \times \vec{A}). \ \ \ \ \ (4)$

Now we discuss under which condition the tight-binding limit can be resonable. The design of the cavities must fulfilled the following three conditions. {}{}

• (1) The resonant frequency of the caivities must falls within the forbidden gap of the surrounding struction. The coupling between cavities is due to the evanescent Bloch waves.
• (2) Cavities are weakly coupled as the cavities are sufficently separated.
• (3) The eigenmode of the field in a cavity of the system will remain nearly the same as the mode in a single cavity. Therefore, only nearest neighbor coupling between cavity modes must be taken into account.

Using tight-binding approximation, ${\vec{A}}$ can be expanded in Wannier function ${\omega_\Omega (\vec{r},t)}$, each located at one single cavity at location ${\vec{R}}$. Besides, it can also be expanded in Bloch wave ${\vec{A}_{\vec{k}} (\vec{r},t)}$. The relation between Wannier function and Bloch function is

$\displaystyle \vec{A}_{\vec{k}} (\vec{r},t) = A_0 e^{i\omega_{\vec{k}} t} \sum_n e^{-ink\vec{R}} \vec{\omega}_{\Omega} (\vec{r} - nR\vec{e}_z). \ \ \ \ \ (5)$

From Eq. (3), we get

\displaystyle \begin{aligned} \nabla \times (\nabla \times \vec{A}_{\vec{k}} (\vec{r},t) ) &= \epsilon( \vec{r}) \frac{\omega_{\vec{k} }^2}{c^2} \vec{A}_{\vec{k}} (\vec{r},t), \\ \nabla \times (\nabla \times \vec{\omega}_\Omega (\vec{r},t) )&= \epsilon_0 (\vec{r},t) \frac{\omega_c^2}{c^2} \vec{\omega}_\Omega (\vec{r},t). \end{aligned} \ \ \ \ \ (6)

${\omega_{\Omega} (\vec{r},t)}$ is real and normalized, ${ \int \mathrm{d}^3 \vec{r} \epsilon_0 (\vec{r}) \vec{\omega}_{\Omega} (\vec{r}) \cdot \vec{\omega}_{\Omega} (\vec{r}) = 1}$. Multiplying both sides of Eq. (6) from left-hand side by Eq. (5) and spatailly integrating, we get the left hand-side of Eq. (6) becomes

\displaystyle \begin{aligned} &\int \mathrm{d}^3 \vec{r} \epsilon(\vec{r}) \frac{ \omega_{ \vec{k}}^2}{c^2} A_0 e^{i\omega_{\vec{k}}t} \sum_n e^{in K \vec{R}} \vec{\omega}_\Omega (\vec{r}) \cdot \vec{\omega}_\Omega (\vec{r} - nR\vec{e}_z) \\ =&\big[ 1 + \Delta \alpha \sum_{n\neq 0} e^{-inkR} \alpha_n \big] \frac{\omega_{k}^2 }{c^2} A_0 e^{i\omega_k t}, \end{aligned} \ \ \ \ \ (7)

where ${\alpha_n = \int \mathrm{d}^3 \vec{r} \epsilon(\vec{r}) \vec{\omega}_{\Omega} (\vec{r}) \cdot \vec{\omega}_\Omega (\vec{r} - n R\vec{e}_z)}$ and ${\Delta \alpha = \int \mathrm{d}^3 \vec{r} [\epsilon(\vec{r}) - \epsilon_0 (\vec{r}) ] \vec{\omega}_\Omega (\vec{r}) \cdot \vec{\omega}_\Omega (\vec{r})}$. The right hand-side of Eq. (6) reads

\displaystyle \begin{aligned} &\int \mathrm{d}^3 \vec{r} \vec{\omega}_{\Omega} (\vec{r}) \cdot \nabla \times (\nabla \times \vec{A}_{\vec{k}} )\\ =&\big[ 1+ \sum_{n\neq 0} e^{-in kR} \beta_n \big] A_0 e^{i\omega_{\vec{k}} t} \frac{\omega_c^2}{c^2} , \end{aligned} \ \ \ \ \ (8)

where ${\beta_n = \int \mathrm{d}^3 \vec{r} \epsilon_0 (\vec{r}- nR \vec{e}_z )\vec{\omega}_\Omega (\vec{r}) \cdot \vec{\omega}_\Omega (\vec{r}- nR \vec{e}_z)}$. Combining Eq. (7) and (8), we get dispersion relation

$\displaystyle \omega_{\vec{k}} \big[ 1+ \Delta \alpha + \sum_{n\neq 0} e^{-in kR} \alpha_n \big] = \omega_c^2 \big[ 1+ \sum_{n\neq 0} e^{-inkR} \beta_n \big]. \ \ \ \ \ (9)$

If the coupling between the resonantor is sufficiently weak, we can count only the nearest neighbor coupling, ie. ${\alpha_n = \beta_n = 0}$, if ${n\neq \pm 1}$. For symmetry considerations, we also require that ${\alpha_1 = \alpha_{-1}}$, ${ \beta_1 = \beta_{-1}}$. Finally we assume ${\alpha_1}$, ${\beta_1}$ and ${\Delta \alpha}$ to be small, we can simplify the dispersion relation to

$\displaystyle \omega_{k} = \omega_c \big[ 1 - \frac{\Delta \alpha}{2} + \kappa_1 \cos (kR) \big], \ \ \ \ \ (10)$

where we define the coupling factor ${\kappa_1}$ as ${\kappa_1 = \beta_1 - \alpha_1}$. Suppose that ${\Delta \alpha}$ is neligiable compared with ${\kappa_1}$. The dispersion relation can be simplified as

$\displaystyle \omega_k = \omega_c [ 1+ \kappa_1 \cos (kR) ]. \ \ \ \ \ (11)$

Now following the usual way to quantize the system from the eigenmodes ${\omega_k}$, we can get Hamiltonian of the system

$\displaystyle H = \omega_c \sum_k b^\dagger_k b_k [ 1+ \kappa_1 \cos (kR)]. \ \ \ \ \ (12)$

From now on, we use natural unit that ${\hbar = 1}$. ${b^\dagger_k}$ and ${b_k}$ is the creation and annihilation operators of a photon occupying the extended eigenmode. It would be conviend to transform the operators in Wannier representation, where the operators describe the localized eigenmodes (cavity modes). The transformation connections is

$\displaystyle b_k = \frac{1}{\sqrt{N}} \sum_n e^{-inkR} a_n, \ \ \ \ \ (13)$

where ${a_n= a(nR)}$. In the localized modes, the Hamiltonian (12) reads

\displaystyle \begin{aligned} H= &\frac{\omega_c}{N} \sum_k \sum_n \sum_{n'} e^{ikR(n-n')} a^\dagger_n a_{n'} \big[ 1 + \kappa_1 \cos(kR)\big]\\ =& \omega_c \sum_{n=1}^N a_n^\dagger a_n + \frac{\kappa_1 \omega_c }{2N} \sum_{n,n'} \sum_{k} e^{ikR(n-n')} a^\dagger_n a_n' (e^{ikR}+ e^{-ikR}) \\ =& \omega_c \sum_{n=1}^N a_n^\dagger a_n + \frac{\kappa_1 \omega_c}{2N} \sum_{n,n'} \sum_k \Big[ e^{ikR(n-n'+1)} a^\dagger_n a_n' + e^{ikR(n-n'-1)} a^\dagger_n a_n' \Big] \\ =& \omega_c \sum_{n=1}^N a_n^\dagger a_n + \frac{\kappa_1 \omega_c}{2} \sum_{n=1}^N [a^\dagger_n a_{n+1} + \mathrm{H.C.} ] + \frac{N\kappa_1 \omega_c}{2} \end{aligned} \ \ \ \ \ (14)

Neglect the total energy term ${\frac{N\kappa_1 \omega_c}{2}}$, we get the effective Hamiltonian is

$\displaystyle H_{eff} = \omega_c \sum_{n=1}^N a_n^\dagger a_n + \frac{\kappa_1 \omega_c}{2} \sum_{n=1}^N [a^\dagger_n a_{n+1} + \mathrm{H.C.} ], \ \ \ \ \ (15)$

where we use the donate ${a_{N+1} = a_N}$.

3. Site repulsion

To generate a repulsion between polaritons that are located in the same cavity, we fill the cavity with 4 level atoms.

In a rotating frame with ${H_0 = \omega_c (a^\dagger a + 1/2) + \sum_{j=1}^N (\omega_c \sigma_{22}^j + \omega_{c} \sigma_{33}^j + 2 \omega_c \sigma_{44}^j)}$, the Hamiltonian of the atoms in the cavity reads

$\displaystyle H_I = \sum_{j=1}^N \big[\varepsilon \sigma_{22}^j + \delta \sigma_{33}^j + (\Delta + \varepsilon) \sigma_{44}^j \big] + \sum_{j=1}^N ( \Omega_L \omega_{23}^j + g_{13} \sigma_{13}^j a^\dagger + g_{24} \sigma_{24}^j a^\dagger + \mathrm{H.C.} ), \ \ \ \ \ (16)$

where ${\sigma_{kl}^j = |k_j\rangle \langle l_j|}$. ${\omega_c}$ is the frequency of the cavity mode. ${\Omega_L}$ is the Rabi frequency of the driving field. ${g_{13},g_{24}}$ are the parameters of the dipole coupling of the cavity mode.

The Hamiltonian truncated to the subspace at most two excitations can be diagonalized. Let’s define the following operators ${ S^\dagger_{12} = \frac{1}{\sqrt{N}} \sum_{j=1}^N \sigma_{21}^j}$ and ${ S^\dagger_{13} = \frac{1}{\sqrt{N}} \sum_{j=1}^N \sigma_{31}^j}$. In the limit ${N\gg 1}$ and excitation is restricted to ${2}$, we get latex \displaystyle \begin{aligned} \big[ S_{12}, S^\dagger_{12} \big] =& \frac{1}{N} \sum_{j=1}^N |1\rangle \langle 1|_j \approx 1 ,\\ \big[ S_{13}, S^\dagger_{13}\big] = & \frac{1}{N} \sum_{j=1}^N |1\rangle \langle 1|_j \approx 1 ,\\ \big[ S_{12}, S^\dagger_{12} \big] = & \big[ S_{13}, S^\dagger_{13} \big] =0 \\ \big[S_{12}^\dagger, S_{13} \big] = &\big[S_{13}^\dagger, S_{13} \big] = 0. \end{aligned} Therefore, ${S_{12}}$ and ${S_{13}}$ fulfill boson commutation and represent two different bosons.

Neglecting two photon detuning and coupling to level ${4}$. Setting ${g_{24}=0}$ and ${\varepsilon =0}$, level ${4}$ decouples from the rest of the two exciation manifold. We assume the number of atoms is large, ${N \gg 1}$. The Hamiltonian (16) in terms of ${S_{12}}$ and ${S_{13}}$ reads

$\displaystyle H_I = \delta S^\dagger_{13} S_{13} + \Omega_L S^\dagger_{12} S_{13} + g S_{13} a^{\dagger} + \Omega_L S_{13}^\dagger S_{12} + g S_{13}^\dagger a, \ \ \ \ \ (17)$

where ${g= \sqrt{N} g_{13}}$. Let’s define operator ${Q_0^\dagger = \frac{1}{B} ( \Omega_L S_{12}^\dagger + ga^\dagger )}$, where ${B= \sqrt{\Omega_L + g^2}}$. ${Q_0}$ is boson which fulfills ${[Q_0, Q^\dagger_0] = 1}$. We can define another boson operator ${P_0^\dagger= \frac{1}{B} ( g S_{12}^\dagger - \Omega_L a)}$ which fulfils ${[P_0,P^\dagger_0]=1}$, ${[P_0, Q_0] = [P_0, Q_0]=0}$. So ${P_0}$ and ${Q_0}$ are two different bosons.

In terms of ${Q_0}$ and ${S_13}$, Hamiltonian (16) reads

$\displaystyle H_I = \delta S_{13}^\dagger S_{13} + B (Q^\dagger_0 S_{13}+ S^\dagger_{13} Q_0). \ \ \ \ \ (18)$

We want use normal transformation to diagonlize ${H_I}$. Let’s define

\displaystyle \begin{aligned} P^\dagger_+ = &\mu S^\dagger_{13} + \nu Q^\dagger_0 \\ P^\dagger_+ =& -\nu S^\dagger_{13} + \mu Q^\dagger_0, \end{aligned} \ \ \ \ \ (19)

where ${\nu}$ and ${\mu}$ are parameters which fulfill ${|\nu|^2 + |\mu|^2 = 1}$. Is is easy to prove that ${[P_+ ,P^\dagger_+] = [P_-, P^\dagger_-] = 1}$ and ${[P_+, P_-^\dagger] = [P_-, p_+^\dagger] =0}$. So ${P_-}$ and ${P_+}$ are two different bosons. ${H_I}$ in terms of ${P^\dagger_+}$ and ${P^\dagger_-}$ reads

\displaystyle \begin{aligned} H_I = &\delta (\mu P^\dagger_+ - \nu P^\dagger_-) (\mu P_+ - \nu P_-) + B\big[ (\nu P^\dagger_+ + \mu P^\dagger_-) ( \mu P_+ - \nu P_-) \\&+ (\nu P^\dagger_+ - \nu P^\dagger_-) ( \nu P_+ + \mu P_-) \big] \end{aligned} \ \ \ \ \ (20)

To diagonlize ${H_I}$, ${\nu}$ and ${\mu}$ must fulfill

$\displaystyle B(\nu^2 - \mu^2) + \delta \mu\nu = 0.$

Solve equation we get ${\frac{\nu}{\mu} = \frac{\delta \pm A}{2B}}$, where ${A = \sqrt{\delta^2 + 4B^2}}$. The diagonlized Hamiltonian is

$\displaystyle H = \mu_0 P^\dagger_0 P_0 + \mu_+ P^\dagger_+ P_+ + \mu_- P^\dagger_- P_-, \ \ \ \ \ (21)$

where ${\mu_0=0}$, ${\mu_+ = (\delta+A)/2}$ and ${\mu_- = (\delta-A)/2}$.

To write the full Hamiltonian ${H_I}$ in the polariton basis, we express operators ${\sum_{j=1}^N \sigma_{22}^j}$ and ${a^\dagger \sum_{j=1}^N \sigma_{24}^j}$ in terms ${P^\dagger_0}$, ${P_+^\dagger}$ and

${P^\dagger_-}$. Let’s suppose that ${|g_{24}|,|\varepsilon|, |\Delta| \ll |\mu_+ |, |\mu_-|}$. Therefore rotating wave approximation is used here. The interaction between ${P_0}$ and ${P_+}$, ${P_-}$ can be negelected. Therefore we have

\displaystyle \begin{aligned} \varepsilon \sum_{j=1}^N \sigma_{22}^j = &\varepsilon S_{12}^\dagger S_{12} \\ \approx & \varepsilon \frac{g^2}{B^2} P^\dagger P_0 \\ g_{24} \sum_{j=1}^N (\sigma_{42} a + \mathrm{H.C.} =& g_{24} (S^\dagger_{14} S_{12} a + \mathrm{H.C.} \\ \approx & - g_{24} \frac{g \Omega_L}{B^2} (S^\dagger_{14} P_0^2 + \mathrm{H.C.}). \end{aligned} \ \ \ \ \ (22)

At last, level ${4}$ can be adiabatically eliminated. The Hamiltonian we get at last is

$\displaystyle H = \kappa (P_0^\dagger)^2 P^2, \ \ \ \ \ (23)$

where

$\displaystyle \kappa = \frac{g_{24}^2}{\Delta} \frac{N g_{13} \Omega_L}{(N g_{13}^2 + \Omega_L^2 )^2}$

. For simplicity, we set ${\varepsilon =0}$.

4. The total Hamiltonian

If we use polariton operators to write the hopping terms in Eq. (1),

$\displaystyle a^\dagger_n a_{n+1} \approx \frac{\Omega_L^2}{B^2} P^\dagger_n P_{n+1}.$

Contribution of different polaritons decouple due to the separation of their energy. There for we get the Hamiltonian ${P_n^\dagger}$ takes the on the form (1), with

$\displaystyle J= \frac{ \omega_C \Omega^2_L}{2(N g_{13}^2 + \Omega_L^2)} \kappa_1.$