issue with basic Hamiltonian with complex values

[josselin....@gmail.com]

Hi everyone

I'm having an issue with a very simple, basic case:

2-level system with constant Rabi frequency; only twist is I'm adding an additional constant phase factor e^(i*phi)to the Rabi frequency.

when the phase factor is 1 (or -1) , I get Rabi oscillation as expected. Obvious.

for an arbitrary phase factor, however, I get something like a 2-level with a constant drive and a decay (spontaneous emission)

I expect the 2 coefficient of |0> and |1> to just pick up this phase difference between them and this constant phase factor to just get change nothing to the probability |<1|psi(t)|0>|^2, in any case not to have this effect.

I even upgraded to use the QobjEvo class to define my Hamiltonian and tried to follow the doc for that.

I'm probably missing something obvious, either in my math or code.

I attached an example notebook. The dict defines the constant phase, the amplitude of the Rabi freq and where the constant drive stops on the interval from 0 to T (total propagation time).

So again, I'm expecting the phase to change nothing in this notebook.

I'm really confused by this ultrasimple thing.

any help much appreciated...

cheers

Josselin

[josselin....@gmail.com]

I'm just pasting the code from the attached nb, in case that's more the habit here.

there is a function for both array-type and function-type Hamiltonian.

%load_ext watermark
import numpy as np
import matplotlib.pyplot as plt
import qutip
%watermark -v --iversion

N = 300
T=1.5
tlist = np.linspace(0, T, N)
ket0 = qutip.basis(2,0)
ket1 = qutip.basis(2,1)

proj10 = qutip.ket2dm(ket0)
proj01 = qutip.ket2dm(ket1)

def Omega_Cst(times, args):
    A = args['ampl']
    fi = args['phi']
    last = args['Nf']
    res = np.empty(len(times), dtype=np.complex128)
    for k in range(last):
        res[k]=A*np.exp(1.j*fi)
    for k in range(last, len(times)):
        res[k]=0
    return res

def Omega_Cst_f(t, args):
    A = args['ampl']
    fi = args['phi']
    last = args['Nf']
    if t<=T*float(last)/N:
        return A*np.exp(1.j*fi)
    else:
        return 0

dict_cst = {'phi': np.pi*0.5, 'ampl':2*np.pi/T, 'Nf':int(N*1.)}

fig, ax = plt.subplots()
ax.set_ylim(0, dict_cst['ampl']*1.05)
om_cst = Omega_Cst(tlist, dict_cst)
ax.plot(tlist, np.absolute(om_cst))

def solve_cst_func(detuning):
    H = qutip.QobjEvo([detuning*qutip.sigmaz(), [qutip.sigmax(), Omega_Cst_f]], args=dict_cst)
    result = qutip.sesolve(H, ket0, tlist, e_ops=[proj10, proj01])
    return result

def solve_cst_arr(detuning):
    H = qutip.QobjEvo([detuning*qutip.sigmaz(), [qutip.sigmax(), Omega_Cst(tlist,dict_cst)]], tlist=tlist)
    result = qutip.sesolve(H, ket0, tlist, e_ops=[proj10, proj01])
    return result

Cst_omega_slv = solve_cst_func(0)

fig,ax=plt.subplots()
ax.set_ylim(0,1)
ax.plot(tlist, Cst_omega_slv.expect[1])

[Simon Cross]

Hi Josselin,

Looking at your code I think you have made a mistake in how you interpret where / how the phase difference can be applied. For the oscillation you add a term A * sigmax * e^(i * phi), but this term is not even guaranteed to be Hermitean.

For example,

q = qutip.sigmax() * np.exp(np.pi * 0.5 * 1j)
print(q)
print("--")
print(q.eigenenergies())
print("--")
print(q.isherm)

You can rotate the phase of the initial ket, but if one applies an operation U to a ket, the equivalent transformation for an operator such as H is U * H * U.dag().

Hope this helps untangle things.

Regards,

Simon

[josselin....@gmail.com]

oh right ! I got it.

thanks a lot Simon! sorry for the beginner's mistake

I knew I was somehow mathematically in the wrong since all methods where failing...

I guess I should have just written on paper what I was doing and that would have shown me !

very nice of you to take the time

regards

Josselin

[Simon Cross]

Pleasure! Glad it helped. :D