How to manipulate the phase of the electromagnetic radiation in the Hamiltonian

Question

The papers describing the Honeywell trapped ion quantum computer say that rotations of the qubit about an axis in the xy plane of the Bloch sphere are accomplished by applying a slightly detuned pair of lasers whose "beat" frequency matches the frequency of the qubit and whose relative phase determines the azimuthal angle of the rotation axis in the sphere. Using the quantum simulation package QuTiP, I set up the Hamiltonian and solved Schrodinger's equation for this situation. I found that when the photon annihilation operators for the two modes have no relative phase between them, the state evolves in such a way as to reflect rotation about the x axis. This is as expected.

However, in order to get a rotation about the y axis, I need to shift the relative phase of the "beat note" and I don't know what the Hamiltonian would look like.

Here is my code for the Hamiltonian without the phase shift:

wa0=0 #Frequencies
wa1=4*np.pi
w_light = wa1*3000
detuning = wa1-wa0
g_light = 0.02*2*np.pi #Coupling
N = 15                 # number of cavity fock states


def hamiltonian(N, g_light, wa0, wa1, w_light, detuning):
    a_light1  = -tensor(destroy(N), qeye(N), qeye(2)) #Annihilation operator mode 1
    a_light2  = tensor(qeye(N), destroy(N), qeye(2))  #Annihilation operator mode 2
    w_light2 = w_light + detuning                     #Frequency mode 2
    sm10 = tensor(qeye(N), qeye(N), Qobj([[0,1],[0,0]])) #Transition from first excited state to ground state
    
    # Hamiltonian without light-matter interaction
    H = w_light * a_light1.dag() * a_light1 + w_light2 * a_light2.dag()*a_light2 + \
        tensor(qeye(N), qeye(N), Qobj([[wa0,0],[0,wa1]]))
    
    # Light-atom interaction
    Hi = g_light*\
        (\
         a_light1.dag() * sm10 + \
         sm10.dag()*a_light1 + \
         a_light2.dag()* sm10 + \
         sm10.dag()* a_light2 \
        )

    #Transform to interaction picture
    h=[[np.exp(1j*wa0),0],[0,np.exp(1j*wa1)]]
    EH=tensor(qeye(N),qeye(N),Qobj(h))
    H=EH*Hi*EH.dag()

    # It has to be Hermitian
    assert H.isherm
    return H

Here is my initial state:

atom_state = fock(2,0)
psi0 = tensor(coherent(N,2), coherent(N,2), atom_states[state])

"Coherent" is, of course, a coherent optical state (eigenstate of annihilation operator)

Answer 1

Simulating a 2 photon Raman transition process is explained very clearly in lectures 9 and 10 here:

http://people.ee.duke.edu/~jungsang/ECE590_01/

Specifically, take a look here: http://people.ee.duke.edu/~jungsang/ECE590_01/ThreeLevelSystems.pdf

and at programming exercise 7.

To summarize the differences between what you have and the required simulation:

1. The driving fields should be treated as classical driving fields and not a quantum mechanical Bosonic mode.
2. The simulation needs to include a finite linewidth for the state
3. It also needs to include detuning, because we are virtually populating the $|2\rangle$ state.

Once you put all these in as described in the lecture notes, you should be able to see the phase dependet Rabi tranisition.

BTW, to get you started, I am including here some basic qutip code that should contain the necessary ingredients. I didn't get the correct Rabi flipping behavior and it's a bit hard for me to tweak it due to lack of time, but I believe it should put you on the right path.

import qutip as qp
import numpy as np

# parameters
w01 = 1
w12 = 4
g02 = 0.01
g12 = 0.01
gamma01 = 0.01
gamma12 = 0.01


phi = 0  # angle between drives

# operators
a02 = qp.Qobj([[0, 0, 1], [0, 0, 0], [0, 0, 0]])
a12 = qp.Qobj([[0, 0, 0], [0, 0, 1], [0, 0, 0]])
a01 = qp.Qobj([[0, 1, 0], [0, 0, 0], [0, 0, 0]])
H0 = qp.Qobj(np.diag([0, w01, w01 + w12]))
psi0 = qp.fock_dm(3, 0)


def drive02p(t, args):
    return g02 * np.exp(-1j*((w01 + w12)*t))


def drive02m(t, args):
    return g02 * np.exp(1j*((w01 + w12)*t))


def drive12p(t, args):
    return g12 * np.exp(-1j*(w12*t + phi))


def drive12m(t, args):
    return g12 * np.exp(1j*(w12*t + phi))


H = [H0,
    [a02, drive02m], [a02.dag(), drive02p],
    [a12, drive12m], [a12.dag(), drive12p]
    ]

# solve
res = qp.mesolve(H,
                 psi0,
                 np.linspace(0, 200, 20000),
                 e_ops=[a01 + a01.dag(), 1j * a01 + -1j * a01.dag(), H0],
                 c_ops=[np.sqrt(gamma01) * a01, np.sqrt(gamma12) * a12]
                )

plt.plot(res.times, res.expect[0], label='X')
plt.plot(res.times, res.expect[1], label='Y')
plt.plot(res.times, res.expect[2], label='Z')
plt.legend()
```

Answer 2

The given answer did not solve my problem. To simulate phase, it was necessary to initialize with a coherent superposition of photon Fock states (number states). These states are eigenvectors of the lowering operator and the way to encode the phase of the simulated light is to choose and eigenvalue with that complex phase. Because the annihilation operator is not Hermitian, its eigenvalues need not be real.

Once this step was taken, it was a simple matter to simulate gate operations that perform a rotation in the Bloch sphere. The quantum simulation library qutip provides a function "coherent" that returns such a state.

The point is that the phase does not go into the Hamiltonian, it goes into the initial state.