Cavity QED with tunneling: an application of adiabatic elimination in quantum trajectory theory
Charles Baldwin Perry Rice





Macklin Quantum Information Sciences
Department of Physics, Miami University, Oxford OH 45056






We model the evolution of a single atom in a cavity interacting with two lasers, one far off resonance which creates an optical potential lattice and one near resonance, which can interact with the atom.  The atom may tunnel between sites in the lattice.  We find the photon counting statistic g^(2)(τ) depends on the tunneling rate.  Additionally we derive methods to work with a Bose-Einstein Condensate in the same situation, in ther Josephson junction regime.  It is shown that oscillations in g^(2)(τ) are directly proportional to the tunneling coefficient. We work in the bad-cavity limit and give an example of the use of adiabatic elimination in quantum trajectory theory


Cavity QED with tunneling: an application of adiabatic elimination in quantum trajectory theory
    Perry Rice
    December 30, 2023
===============================================================================================








§ INTRODUCTION


Cavity quantum electrodynamics (cQED) is the study of ultra-cold atoms interacting with laser light in a high quality mirror cavity.  


In this paper, we model a single atom in an optical lattice, within a cavity, allowed to tunnel between lattice sites.  The model provides a base for further work with multiple atoms, which will be able to interact.  We treat the atoms as a Bose-Einstein condensate which will be discussed later. We establish methods to study the particle's behavior by observing emission from the cavity.


There have been many variations of the cQED system studied.  We will explain the basic model here.  We place the atoms inside a high finesse mirror cavity.  The atoms are assumed to only have two allowed energy levels, excited (|e⟩) and ground (|g⟩).  The cavity is then pumped with a single mode driving laser of strength Y that interacts with the atoms with coupling rate g.  The atoms will then oscillate between its two energy levels with the Rabi frequency equal to g √(N) for N atoms.

This alone, however, is not a full description of the system.  We have failed to include the loss of photons from the cavity.  There are two types of photon loss, fluorescence and transmission.  Fluorescence is a result of the atom absorbing a photon from the driving laser and later emitting it out the side of the cavity in a random direction (spontaneous emission.)  
 In addition to these interactions we add, in this paper, add an optical lattice and tunneling.  The lattice is created by pumping a far off resonant laser into the cavity creating a finite number of sites that the atom can occupy.  The atoms can jump between wells determined by a constant T.  We can also treat the atoms as  Bose-Einstein condensate, which requires accounting for additional interactions between atoms.







To solve for the evolution of the system we utilize quantum trajectory theory. A non-Hermitian Hamiltonian accounts for all of these interactions modeling the atoms evolution in the cavity.  Quantum jumps occur which are random events modeling the fluorescence and transmission discussed earlier.  Using times when photons are fluoresced (spontaneous emission) and transmitted (cavity loss.)  we can calculate the photon counting statistics.




In section two, we will outline the system in question.  In section three, we will derive the weak field trajectory model which treats atoms in a cavity interacting with a driving field  and an off resonance lattice beam.  In section four, we discuss the manner in which one uses adiabatic elimnation of variables.  I In section five we will combine the weak field trajectory model and the Bose-Einstein condensate model to give a single Hamiltonian to model N two level boson atoms in an optical lattice interacting with a driving field.  In section six, we will present the numerical model of a single two level atom in a two site optical lattice and compare it to theoretical predictions.  And finally, in section seven we conclude the analysis and give ideas for future studies.





§  OPTICAL LATTICES
 
Eventually our system will include a 1D optical lattice that will restrict the motion of the atoms inside the cavity.  Optical lattices represent a periodic potential created by a set of interfering laser beams.  The potential is caused by the ac-stark shift on neutral atoms <cit.>.  The result is the atoms are evenly distributed like a crystal  <cit.>.  However, unlike crystals they are easy to manipulate since the size and strength of the lattice sites are determined by the laser creating the lattice. 

The model we use here is semi-classical where we treat atom quantum mechanically but the field classically.  The Hamiltonian then contains two terms, one for the atom and the other for the atom's dipole interaction with the field

    Ĥ = ħω_0/2σ̂_z - E⃗(z, t) ·d̂

The electric field from a laser is E⃗(z, t) = E_0(z) cos(ω t) x̂ oscillates in x̂ and polarized in x̂, which makes the Hamiltonian

    Ĥ = ħω_0/2σ̂_z - ħΩ(z) cos(ω t) σ̂_x

where Ω = E_0(z) ⟨e | x̂·d̂ | g|$⟩ andσ̂_x = σ̂_+ + σ̂_-which results from the dot product. The Hamiltonian can be rewritten with the rotating wave approximation, which neglects rapidly oscillating terms

    Ĥ≈ħω_0/2σ̂_z - ħΩ(z)/2( e^-i ω tσ̂_+ + e^i ω tσ̂_- )
eq: 2comp BEC Ham
From here we make the transformation to a rotating frame at rate of the laserω.

    Ĥ≈  -ħΔ/2σ̂_z - ħΩ(z)/2σ̂_x

WhereΔ= ω- ω_0is the detuning.  The eigenvalues for the Hamiltonian are the energy levels of the system

    E_± = ±√((ħΔ/2)^2  + ( ħΩ(z)/2)^2)≈±( ħΔ/2 + ħΩ^2(z)/4 Δ)

This is the energy of a two level atom either in the exited (+) or ground state (-).  The first term is constant so we can shift the zero point energy to that value.  The second term varies with position.  Since the lattice is formed by a laser we assume it has the formΩ(z) = Ω_0 sin(kz)which makes the second term proportional to

    E_ls∝ sin^2(kz)

where the subscript denotes that this is the light shift energy <cit.>.  The light shift is proportional to the potential the atom feels in the optical lattice.  We stated before that the lattice is periodic and this derivation verifies that claim.  The laser's wavelength and intensity define the site spacing and depth respectively.  For our model, we treat the lattice like an external potential.

    V(z) = V_0 ^2(kz)

All of the atoms within the cavity feel an external force dependent on their position.  Later we will show how this external potential affects the behavior of a BEC.





§  STATE AMPLITUDES AND WEAK FIELD LIMIT

In this section, we will combine all of our analysis thus far to find the state amplitudes of the system.  We will put a limit to the number of photons in the system, known as the weak field limit, in order to keep the basis a manageable size.  The analysis here follows <cit.> and <cit.>.  

    Ĥ =  i ħΔ_a/2Ĵ_z + ħΔ_c â^†â +  i ħ g (  â^†Ĵ_- - âĴ_+ ) + i ħ Y ( â^† -  â)- iħ (γ/2Ĵ_+ Ĵ_- + κâ^†â)

HereΔ_a = ω- ω_0andΔ_c = ω- ω_care the detuning between the atom's internal frequency (ω_a) and the cavity mode (ω_c) and the driving frequency (ω) respectfully <cit.>.   We add the fourth term to represent a driving laser off resonance.  Theσ̂have also been replaced withĴoperators which are generalized for multiple atoms.  

    Ĵ_± = ∑_i = 1^N σ̂_±^i

UsingĴallows us to look at systems withNatoms which is part of our ultimate goal since dealing with a BEC requires having multiple atoms.  They act on the symmetric states


    |1⟩_S    =   1/√(N)∑^N_k=1|↓_1 ↓_2 ⋯↓_k-1↑_k ↓_k+1⋯↓_N⟩
    |2⟩_S    =   2/√(N(N-1))∑^N_k=1∑^N_l=1, l ≠ k|↓_1 ↓_2 ⋯↓_k-1↑_k ↓_k+1⋯↓_l-1↑l↓l+1⋯↓_N⟩


These states are necessary since we are not able to determine specifically which atom is in the ground state or the excited state <cit.>.  Applying theĴ_±to the symmetric states gives the following
 

    Ĵ_+ |0⟩_S    =   |1⟩_S 
    Ĵ_- |1⟩_S    =   |0⟩_S


And the same can be done with|2⟩_SThe Hamiltonian we will deal with from now on is

    Ĥ =   i ħ g (  â^†Ĵ_- - âĴ_+ ) + i ħ Y ( â^† -  â)- iħ (γ/2Ĵ_+ Ĵ_- + κâ^†â)

If we were to solve the Schroedinger equation from this Hamiltonian there would be an infinite number of states.  In order to truncate the basis we apply the weak field limit where there is less than two excitations (the number of photons plus the number of excited atoms is equal to the number of excitations.)
  

    Ċ_00   =    -Y C_10
    Ċ_10   =    Y (C_00 - C_20)+ g√(N) C_11 - κ C_10
    Ċ_01   =    Y C_11 - g√(N)  C_00 - γ/2 C_01
    Ċ_20   =   √(2) Y C_10 +g √(2)√(N) C_11 - 2 κ C_20
    Ċ_11   =     Y C_01 - g √(2)√(N) C_01 + g √(2)√(N-1) C_02 - (κ + γ/2) C_20
    Ċ_02   =    - g √(2)√(N-1)C_11 - γ C_02

  
HereC_nmare the amplitudes of the|n m⟩state wherenis the number of photons andmis the number of atoms in the excited state.  In the weak field limit, we haveC_nm≈Y^(n+m)/2. The excitation of the system is then referred to as the sumn + m.

This is, however, still not correct since we are missing higher order terms, with the truncated basis, that should appear withĊ_20,Ċ_11, andĊ_02.  In order to fix this problem we go back to the initial assumption of a weak field limit.    Therefore some terms multiplied byY^s, wheresis some integer dependent on the state, do not affect the amplitude rates.  Then we can scale the amplitudes to be on order 1,Y, orY^2grouped by the excitation <cit.>.  We assume the ground stateC_00 ∼1therefore we drop theC_10term inĊ_00since its at least of orderY.  From here it is clear to see that the first excited states (C_10andC_01) are orderYand then we can drop any higher order terms sinceYis small) <cit.>.  Similarly the second excited states (C_20,C_11andC_02) are of orderY^2and so higher orderYterms can be dropped.  The final states are 
   
    Ċ_00   =    0
    Ċ_10   =    Y + g√(N) C_11 - κ C_10
    Ċ_01   =     - g√(N)  C_00 - γ/2 C_01
    Ċ_20   =   √(2) Y C_10 +g √(2)√(N) C_11 - 2 κ C_20
    Ċ_11   =     Y C_01 - g √(2)√(N) C_01 + g √(2)√(N-1) C_02 - (κ + γ/2) C_20
    Ċ_02   =    - g √(2)√(N-1)C_11 - γ C_02


These are the final amplitudes for the weak field trajectory model <cit.>.  We must remember, though, that in this weak field limit, where the amplitudes are scaled to 1,Y,Y^2,..., that all other quantities calculated from them must also be scaled <cit.>.  For example, the jump probability discussed earlier.  



To make our model more efficient we will need to make some approximations.  First, we will describe adiabatic elimination to further simplify our Hamiltonian.  Next, we introduce photon counting statistics which will be our main tool in investigating the system.  



§ ADIABATIC ELIMINATION

Adiabatic elimination is a method to determine the behavior of the system in different time intervals by invoking limits to the parameters involved.   We start here from the Hamiltonian used in the previous section, equation (<ref>).  The limit we will use in this paper is the bad cavity limit meaningκdominates all other rates.  

We have
   
    â̇   =   i/ħ [ H,â] = - gĴ_- - Y + κâ
    â̇^†   =   i/ħ [ H,â^†] =  gĴ_- + Y - κâ^†


The equations forσ̇̂̇_-andσ̇̂̇_+can be solved similarly, however, that is not the simplification we wish to study here <cit.>.  Taking the limit where the rates ofâ̇ →0andâ̇^† →0we can get the equations forâandâ^†. 


    â   =   g/κĴ_- + Y/κ
    â^†   =   g/κĴ_+ + Y/κ


The Hamiltonian we derived previously comes from the density operator <cit.>.  As explained earlier the system is interacting with an external reservoir, therefore we cannot simply replace terms in equation (<ref>) with our new expressions forâandâ^†.  In order to derive the correct Hamiltonian we must start back with the density operator 

    ρ̇ = -i/ħ[ H, ρ] + γ/2[2 Ĵ_- ρĴ_+ - Ĵ_+ Ĵ_-  ρ - ρĴ_+ Ĵ_- ] + κ[ 2 âρâ^† - â^†âρ -  ρâ^†â]

Now we can substitute the newâandâ^†into the master equation.  Instead of substituting in the terms all at once, we will go through each term that is changed starting first with the interaction term in the Hamiltonian.  

The first part comes from the interaction term

    i ħ g (â^†Ĵ_- - âĴ_+ )    =    i ħ g [ ( g/κĴ_+ + Y/κ) Ĵ_- - (g/κĴ_- + Y/κ) Ĵ_+] 
       =    i ħ g Y/κ (Ĵ_- - Ĵ_+) + i ħg^2/κ (Ĵ_+ Ĵ_- - Ĵ_- Ĵ_+)

In this part the second term is equivalent toĴ_zand can be ignored since it only changes the zero point energy.  

The next term changed is the weak driving laser

    i ħ Y (â^† - â) = i ħ g Y/κ (Ĵ_+ - Ĵ_-) + i ħY^3/κ^2

Again, since the final term is higher, it can be ignored.  The first term is equal but opposite to the first term coming from the interaction term and so, they eliminate each other.  This may be disconcerting at first but we must remember that there are more parts from theκportion of equation forρ̇.  

Theκdissipation term has three parts from the foiling of the equations forâandâ^†κ[ 2 âρâ^† -  â^†âρ - ρâ^†â]    =   g^2/κ(2 Ĵ_- ρĴ_+ - Ĵ_+ Ĵ_-  ρ - ρĴ_+ Ĵ_-)
           + i/ħ[i ħg Y/κ (Ĵ_- - Ĵ_+), ρ]

The first part comes from theg Ĵ/κwhile the second comes from the cross term between theĴand the constant term.  There was a third term, which was constant, which we again eliminate by shifting the zero point energy.  The first term has the same operation onρas theγterm and so we can combine the two to give a new decay term in our Hamiltonian or lettingγ→γ/2( 1 + 2 g^2/κγ).  We move the second part to the evolution Hamiltonian to create the bad cavity Hamiltonian

    Ĥ_BC = i ħg Y/κ(Ĵ_- - Ĵ_+ ) - γ/2( 1 + 2 g^2/κγ) Ĵ_+ Ĵ_-

This Hamiltonian is much simpler to work with since we eliminated the field variables.  By using this method, we can solve analytically for some situations.  The probability of a jump from transmission is also edited because of the substitution forâandâ^†as well as the collapse operator.

    P_κ = 2 κ( Y^2/κ^2 + g Y/κ^2 (⟨Ĵ_+|+⟩⟨Ĵ_-|)⟩+ g^2/κ^2⟨Ĵ_+ Ĵ_-|⟩) Δ t

Each term corresponds to a different type of emission <cit.>.  The first is from photons passing through the cavity without being absorbed by an atom.  The second term is from coherent emission by the atoms into the cavity mode.  The third term is from cavity enhanced spontaneous emission of the atoms into the cavity mode.






§ TWO-COMPONENT BEC MODEL

We can now combine the weak field trajectory terms, equation (<ref>), with the Bose-Hubbard Hamiltonian, to get a new Hamiltonian describing the evolution of a two level boson atom in an optical lattice.  The Hamiltonian for the system is

    H    =    - ħ J ∑_i,jb̂^†_i b̂_j + U/2∑_i n̂_i (n̂_i -1 ) +  i ħ g (  â^†Ĵ_- - âĴ_+) + i ħ Y ( â^† -  â) 
           - iħ(γ/2Ĵ_+ Ĵ_- + κâ^†â)

The Hamiltonian allows boson atoms to travel between lattice sites and also interact with a weak driving laserY.  

The number of states for this system is the product of the number of states from the Bose-Hubbard model and the the number of states from the trajectory  model.  The number of states from the Bose-Hubbard Hamiltonian is proportional to the number of sites,Land the number of particlesN<cit.>
    S_BH = (N + L - 1)!/N!(L - 1)!

The number of states for the trajectory Hamiltonian is proportional to the number of particles and the number of excitations.  In the weak field limit there are six possible states forNatoms and five for a single atom.   Then the total number of states is 

    S_tot=6 ×(N + L - 1)!/N!(L - 1)!

For example, when there are two sites and one particles there are two states from the atom configuration (|10⟩and|01⟩where a 1 denotes the one atom in a site) multiplied by the different trajectory states which for the weak filed limit (two excitations) is five.  This gives ten possible states. 

There is also the bad cavity Hamiltonian found by adiabatic elimination of aâandâ^†of the trajectory model as shown in section (<ref>). 

    H_BC =- ħ J ∑_i,jb̂^†_i b̂_j + U/2∑_i n̂_i (n̂_i -1) +  i ħg Y/κ(Ĵ_- - Ĵ_+ ) - γ/2( 1 + 2 g^2/κγ)Ĵ_+ Ĵ_-

In this case the number of states based on the configuration of atoms is the same as the previous Hamiltonian, however, the number of states from the bad cavity Hamiltonian is only equal to the number of excited states sinceâandâ^†were eliminated.  

    S_tot= (e_max + 1) ×(N + L - 1)!/N!(L - 1)!

Where the "+1" represents the ground state.  For the weak field limite_max =2forNatoms while for a single atome_max=1.   Therefore, for the two sites one particle arrangement there are only four possible states.  The number of states from the atom configuration are multiplied by two now.  For the following sections we will restrict our analysis to this system of a single atom allowed to jump between two lattice sites as shown in figure (<ref>).





Now we analyze the single atom two-cavity system with the bad cavity Hamiltonian.  We use the bad cavity Hamiltonian to solve Schroedinger's equation to get the rate of change of the amplitudes.  TheUterm is dropped since that term only appears when there are more than one atoms.  This also means that we are not considering an actual BEC since there are not multiple atoms to interact. 

We use the notationC_nm^iwhere the subscript shows the number of atoms in site one (n) and the number of atoms in site two (m) and the superscript denotes the internal energyi=g,e.  Solving Schroedinger's equation gives


    Ċ_10^g    =    i J C_01^g + g Y/κ C_10^e  
    Ċ_01^g    =    i J C_10^g + g Y/κ C_01^e 
    Ċ_10^e    =    i J C_01^e - g Y/κ C_10^g - γ/2 (1 + 2 C) C_10^e 
    Ċ_01^e    =    i J C_10^e - g Y/κ C_01^g - γ/2 (1 + 2 C) C_01^e


WhereC = g^2/κγ.  From here we can solve the equations numerically or take a weak field limit and solve them analytically.  The numerical solution requires a Runge-Kutta method instead of an Euler method.

The amplitude rates are scaled onYsimilar to in the simple weak field trajectory model shown in section (<ref>.)  The ground states are of order 1 and the excited states are of orderY.  This results in dropping theYterms in the ground state rates


    Ċ_10^g    =    i J C_01^g   
    Ċ_01^g    =    i J C_10^g 
    Ċ_10^e    =    i J C_01^e - g Y/κ C_10^g - γ/2 (1 + 2 C) C_10^e 
    Ċ_01^e    =    i J C_10^e - g Y/κ C_01^g - γ/2 (1 + 2 C) C_01^e


These equations are coupled so to separate them we use the symmetry states labeledD^g,e_±which relate toC^g,eby


    D_+^g    =   1/√(2) ( C_10^g + C_01^g) 
    
    D_-^g    =   1/√(2) ( C_10^g - C_01^g) 
    
    D_+^e    =   1/√(2) ( C_10^e + C_01^e) 
    
    D_-^e    =   1/√(2) ( C_10^e - C_01^e)


Substituting theDequations into equation (<ref>) gives a series of uncoupled differential equations which can be solved analytically.  Undoing the substations forDgives the amplitudes for the weak field bad cavity limit as


    C_10^g(t)    =   1/√(2) (D_+^g(0) e^i J t + D_-^g(0) e^-i J t) 
    
    C_01^g(t)    =   1/√(2) (D_+^g(0) e^i J t - D_-^g(0) e^-i J t) 
    
    C_10^e(t)    =   ( Y D_+^g(0)/Γ (e^-Γ t -1) + D_+^e(0) e^-Γ t) e^i J t/√(2)
           + ( Y D_-^g(0)/Γ (e^-Γ t-1) + D_-^e(0) e^-Γ t) e^-i J t/√(2)
    
    C_01^e(t)    =   ( Y D_+^g(0)/Γ (e^-Γ t-1) + D_+^e(0) e^Γ t) e^i J t/√(2)
          -  ( Y D_-^g(0)/Γ (e^-Γ t-1) + D_-^e(0) e^-Γ t) e^-i J t/√(2)


We letγ(1 + 2C)/ = Γ.  These amplitudes only apply whenYis small due to the weak field limit.  Therefore, we have two separate models to compare, one numeric and one analytic.  To simplify the solution we can restrict the analysis by applying the initial condition that the particle must start in the ground state of the first well.  Then,D_+^g(0) = D_-^g(0) = 1/√(2)andD_+^e(0) = D_-^e(0) = 0.  Then the amplitudes are


    C_10^g(t)    =   (J t) 
    
    C_01^g(t)    =   (J t) 
    
    C_10^e(t)    =   Y/Γ( e^-  Γ t - 1 ) (Jt) 
    
    C_01^e(t)    =   Y/Γ( e^-  Γ t - 1 ) (J t)


Starting the particle in the other well will only shift the phase of each byπ/2.  This will be the starting condition used throughout.  



We start the system in the ground state with the atom in the first well.  The program then does a series of stpdf at a set time interval.  At each step the probability of a jump is calculated,  for spontaneous emission or cavity loss.  For this case (2 sites 1 atom,) there are three jumps.  Two for spontaneous emission (one for each sites) and one for cavity loss.  We can also look at a jump from either site, which is the sum of the jumps from each site.
  

    P_γ_1   =   γ | C_10^e |^2 Δ t 
    
    P_γ_2   =   γ | C_01^e |^2 Δ t 
    
    P_γ   =    P_γ_1 + P_γ_1 =  γ( | C_10^e |^2 + | C_10^e |^2 ) Δ t 
    
    P_κ   =    2 [ Y^2/κ + Y g/κ (C_10^e^*C_10^g + C_01^e^*C_01^g + C_10^g^*C_10^e + C_01^g^*C_01^e) 
          +g^2/κ ( |C_10^e |^2 + |C_01^e|^2)] Δ t


The program then generates a random number and if it falls into a range determined by the probabilities calculated above a jump occurs.  The time of the jump is stored and the collapse operators collapse the system


    σ̂_-^1 |ψ⟩   =    C_10^e |10g⟩
    σ̂_-^2 |ψ⟩   =    C_01^e |01g⟩
    â|ψ⟩   =   (Y/κ C_10^g + g/κ C_10^e) |10g⟩ +  (Y/κ C_01^g + g/κ C_01^e) |01g⟩
          +Y/κC^e_10|10e⟩ + Y/κ C^e_01|01e⟩


If however the random number does not fall into the jump range then the system will evolve.  The numeric and weak field analytic methods differ in this evolution.  For the numerical evolution, the program uses the Schroedinger equation solutions with a Runge-Kutta method.  The analytical model uses the solutions for the amplitudes.  With either case the final amplitudes are renormalized since the Hamiltonian's are non-Hermitian.  
Ultimately, the data taken will be the times of the jumps but to check the model we can  the expectation values⟨σ̂_+^i σ̂_-^i|$⟩ and ⟨â^†â|$⟩ which correspond to the probability of each jump.


    ⟨σ̂_+^1 σ̂_-^1|   ⟩=    |C_10^e|^2 
    ⟨σ̂_+ σ̂_-|   ⟩=     |C_10^e|^2 +  |C_01^e|^2 
    ⟨â^†â|   ⟩=   Y^2/κ + Y g/κ (C_10^e^*C_10^g + C_01^e^*C_01^g + C_10^g^*C_10^e + C_01^g^*C_01^e) 
          +g^2/κ ( |C_10^e |^2 + |C_01^e|^2)


Which are plotted below withJ=1,g=√(10),κ=10,γ=1soC=1.  Each type of jump is labeled to show how the collapse affects the expectation value of the function.