Dear QuTip enthusiasts,
first of all: Thank you very much for your nice work on the Toolbox. It has been fun to work with it so far, and now I face a problem that might also be both interesting and instructive for the community. For the sake of clarity, I will try to structure it in PARTs. For the interested reader, the newest code is also attached – and will be updated to be available as an example, if we are able to fix the problems. :)
PART1 - The problem:I would like to determine the driving field (Omega) that is required to store a photon wavepacket (Phi) in a three-level atom (g1,g2,e) coupled to a cavity (g). For this, the GRAPE algorithm should be used to find, without any initial guess, the required pulse shape Omega for a given photon pulse shape Phi. An analytical solution has been suggested in https://arxiv.org/abs/1105.1699 and is implemented in my jupyter notebook. It should be possible to reproduce that result by solving the problem “just” numerically using QuTip… :)
However, I’m stuck in the definition of the “drive” and “ctrls” terms for the function “optimize_pulse(drift, ctrls,…)”. Following https://github.com/qutip/qutip-notebooks/blob/master/examples/control-pulseoptim-Lindbladian.ipynb , the inputs are created as shown in PART3. The output of the time-dependent driving field amplitudes are always 0, i.e. no pulse is obtained.
PART2 - The system’s Hamiltonian & Collapse operators:
(The rates are scaled such that approximately g = 1.)
#Construct composite operators:
idatom = qeye(3)
idcavity = qeye(N)
g1, g2, e = qutrit_basis()
a = tensor(idatom, destroy(N))
g1_g1 = tensor(g1*g1.dag(), idcavity)
g2_g2 = tensor(g2*g2.dag(), idcavity)
g1_g2dag = tensor(g1*g2.dag(), idcavity)
g2_g1dag = tensor(g2*g1.dag(), idcavity)
e_e = tensor(e*e.dag(), idcavity)
g1_e = tensor(g1*e.dag(), idcavity)
g2_e = tensor(g2*e.dag(), idcavity)#Write down time-independent Hamiltonian-terms:
H_g = g*(1j*a*g2_e.dag() - 1j*a.dag()*g2_e)# Write down time-dependent Hamiltonian-terms
H_Omega = 0.5*(1j*g1_e.dag() - 1j*g1_e)
H_Phi = (a.dag()+a)### combine the time-independent H's ####
H_no_t = H_g# Collapse operators are:
C_a = sqrt(2 *kappa)*a # cavity decay
C_g1 = sqrt(2 *Gamma1)*g1_e # excited state decays
C_g2 = sqrt(2 *Gamma2)*g2_e
C = [C_a, C_g1, C_g2]# The system may then properly described for an initial state psi0 and the driving terms
#“H_Omega,Phi_coeff” using:
H = [H_no_t, [H_Omega, H_Omega_coeff], [H_Phi, H_Phi_coeff]]
output = mesolve(H, psi0, tlist, C, [a.dag()*a, g1_g1, g2_g2, e_e], args=args)
PART3 – Optimal control
As described in PART1, I followed the example setting:
drift = list(map(lambda t: (liouvillian(H_no_t + H_Phi_coeff(t, args) * H_Phi, [C_a,C_g1,C_g2])), tlist))
ctrls = liouvillian(H_Omega),since drift has to be an array of Hamiltonians, if we regard a time-dependent driving of the cavity. If here also lies a mistake, I would be interested in how to set up the problem at all, since Phi cannot go into the ctrls terms that will be altered. I noticed that moving the Collapse operators also to ctrls, a – not yet optimal - output is created. But why would I have to split or two-fold implement them?
Part4 – Python and Qutip version
umpy Version: 1.14.5
Scipy Version: 1.1.0
Cython Version: 0.28.3
Matplotlib Version: 2.2.2
Python Version: 3.6.5
Number of CPUs: 12
BLAS Info: OPENBLAS
OPENMP Installed: False
INTEL MKL Ext: False
Platform Info: Linux (x86_64)
I’d be very happy to hear your thoughts on this problem and to answer any questions!
All the best,
Tobias