Strange behavior of mesolve for a time depended problem in the presence of collapse operators

[nlphysics]

Hi,

I was writing a simple code for a time dependent driving of a two level atom with the excited state decay and encountered the following problem: If I evolve the system using mesolve for particular time dependence of my driving field, namely if the evolution time is set to be from t = 0 to t = 3.5 in the units of 1/(Rabi frequency) and the driving is only present between t = pi/2 and t = pi, I get expected results if I do not include any collapse operators (using []), however as soon as I use the usual collapse operator corresponding to the excited state decay it gives me for all the time the initial values of the quantities. It happens even if I set the decay rate to zero! Also, mcsolve works fine in both cases. Here is the code that I use:

import numpy as np
import pylab as plt
import scipy as scipy
from qutip import *

#System parameters
Omega = 1
Gamma = 0.1

g = basis(2,0)
e = basis(2,1)
sgg = g * g.dag()   #ground state population
see = e * e.dag()   #excited state population
sge = g * e.dag()   #atomic coherence sigma_ge

#Driving Hamiltonian
H1 = Omega * (sge + sge.dag())

#Time dependence of the drive
def H1_coeff(t,args):
    if t > np.pi:
        c =0
    elif t >= np.pi/2:
        c = 1
    else:
        c = 0
    return c

#Complete Hamiltonian
H = [[H1, H1_coeff]]

#Time
tspan = np.linspace(0,3.5,500)

#initial state
psi0 = g

#collapse operators
c_ops = [np.sqrt(Gamma) * sge]

#calculated expectation values
e_ops = [sgg, see]
results = mesolve(H, psi0, tspan, c_ops, e_ops)

#plotting the results
fig, ax = plt.subplots(1,1)
ax.plot(tspan, results.expect[0])
ax.plot(tspan, results.expect[1])
plt.show()

And here is the information about the packages that I'm using

QuTiP: Quantum Toolbox in Python
Copyright (c) 2011 and later.
A. J. Pitchford, P. D. Nation, R. J. Johansson, A. Grimsmo, and C. Granade

QuTiP Version:      4.2.0
Numpy Version:      1.13.1
Scipy Version:      0.19.1
Cython Version:     0.26
Matplotlib Version: 2.0.2
Python Version:     3.6.2
Number of CPUs:     2
BLAS Info:          INTEL MKL
OPENMP Installed:   False
INTEL MKL Ext:      True
Platform Info:      Darwin (x86_64)

Another funny fact that I noticed is that changing the elif condition to t >= 1.4 gives nice results but changing it to t >= 1.5 will again deliver only the initial values.

[nonhermitian]

My guess is it is a side effect of using the step function time-dependence. Instead, try a tanh, or other continuous and differentiable approximation, and see if that works.

[nlphysics]

That was my first thought as well and it worked for the above example. However, it seems that this issue somehow depends on the dimensionality of the system, since what I actually need is a system of two three level atoms which can be driven independently. So if I modify the code to

g = basis(3,0)
e = basis(3,2)
s = basis(3,1)

#Single atom operators
sgg = g * g.dag()   #ground state population
see = e * e.dag()   #excited state population
sge = g * e.dag()   #atomic coherence sigma_ge

#Composed atomic operators
A1sgg = tensor(sgg, qeye(3))
A1see = tensor(see, qeye(3))
A1sge = tensor(sge, qeye(3))

H1 = Omega1 * (A1sge + A1sge.dag())
def H1_coeff(t,args):
    return 0.5 * (np.tanh(100.0*(t-np.pi/2.0))-np.tanh(100.0*(t-np.pi)))
H = [[H1, H1_coeff]]
c_ops = np.sqrt(Gamma) * A1sge
e_ops = [A1sgg, A1see]
psi0 = tensor(g, g)

I again get only the initial values. With that being said, the same time dependence function works fine for two two level atoms. So it's still puzzling to me. And again if I put [] for the c_ops it works! After looking into mesolve source code I noticed that in that case other solver function is used sesolve, which uses the states without converting them to density matrices and hence has lower dimensionality.

[nonhermitian]

Did you try decreasing Gamma?

[nlphysics]

Yes! I played with Gamma as well as with the number of time steps. Nothing helped. As I already mentioned in the previous post even setting Gamma to zero doesn't help. What I noticed though if I change the factor within the tanh function from 100 to 90 it works again. It's really weird, since it worked with factor 100 in the case of two two level atoms but not for two three level atoms even after increasing the number of time steps by a factor of 20. I really don't know what's going on.

[nonhermitian]

The issue still is the rapid slope on the edges of the pulse. The default tolerance on the ODE solver is not enough to capture the dynamics in your case. Doing:

options=Options(atol=1e-9)

in mesolve resolves the issue. The evolution in the sesolve case is a bit different since the vector being solved is a ket rather than dm.

On Oct 19, 2017, at 22:06, nlphysics notifications@github.com wrote:

Yes! I played with Gamma as well as with the number of time steps. Nothing helped. As I already mentioned in the previous post even setting Gamma to zero doesn't help. What I noticed though if I change the factor within the tanh function from 100 to 90 it works again. It's really weird, since it worked with factor 100 in the case of two two level atoms but not for two three level atoms even after increasing the number of time steps by a factor of 20. I really don't know what's going on.

— You are receiving this because you were assigned. Reply to this email directly, view it on GitHub https://github.com/qutip/qutip/issues/771#issuecomment-338100758, or mute the thread https://github.com/notifications/unsubscribe-auth/ABMPqTRV0BuVMs-vNMMGSkAmYK_qmopBks5suBxXgaJpZM4P__42.

[nlphysics]

After setting atol=1e-9 it works now for 500 time steps. But increasing the number of time steps to 550 causes the program to fail again. It also doesn't help to decrease atol further, I went up to 1e-15, neither does decreasing the relative tolerance rtol. Shouldn't actually increasing the number of time steps improve the precision? Also, why does it depend on the dimension of the dm? I'm just curious now.

[nonhermitian]

There is also an rtol value to play with. However, when the pulse is off, your Hamiltonian is a zero operator. Shouldn't you have the bare terms that give the energies of your states?

[nlphysics]

As mentioned in the previous post I also played with rtol value, but it didn't help. Concerning your question, yes when the pulse is off, I have a zero Hamiltonian and the system stays in the initial state, which is fine, the weird thing however is that apparently the system stays in the initial state even when the pulse is on. Basically, I get as a result that the system stays in the initial state for all times. Did you try to run the code on your machine? Maybe it's the issue of my packages configuration?

[nonhermitian]

I have not run the latest incarnation of your code. However, I am pretty sure it is not a QuTiP issue, but something to do with your implementation. Again, the use of a zero Hamiltonian is suspect. It also doesn't make sense with the rest of your code where you assume ground, excited, and an 's' state, but your Hamiltonian is degenerate. Since you are free to add a constant to any Hamiltonian, I would at least not use a zero operator.

[nlphysics]

Adding a constant Hamiltonian did also not help. In my most recent code I even removed all the atomic states degeneracies and it still won't work. I post my code below. Could you run it on your machine and let me know, whether you get more than constant lines as an output.

import numpy as np
import pylab as plt
import scipy as scipy
from qutip import *

#System parameters
Omega = 1
Gamma = 0.25

#Single atom states
g = basis(3,0)
e = basis(3,2)
s = basis(3,1)

#Single atom operators
sgg = g * g.dag()   #ground state population
see = e * e.dag()   #excited state population
sss = s * s.dag()
sge = g * e.dag()   #atomic coherence sigma_ge
#Composed atomic operators
A1sgg = tensor(sgg, qeye(3))
A1see = tensor(see, qeye(3))
A1sge = tensor(sge, qeye(3))
A1sss = tensor(sss, qeye(3))
A2sgg = tensor(qeye(3), sgg)
A2see = tensor(qeye(3), see)
A2sss = tensor(qeye(3), sss)

#Time independent part of the Hamiltonian
H0 = -1 * A1sgg + 0.5 * A1sss + 1 * A1see - 0.5 * A2sgg + 0.1 * A2sss + 1.2 * A2see
#Driving atom A1
H1 = Omega * (A1sge + A1sge.dag())
#Time dependence of the driving
def H1_coeff(t,args):
    return 0.5 * (np.tanh(100.0*(t-np.pi/2.0))-np.tanh(100.0*(t-np.pi)))
#Complete Hamiltonian
H = [H0, [H1, H1_coeff]]

#Time
tspan = np.linspace(0,3.5,550)
#initial state
psi0 = tensor(g,g)
#calculated expectation values
e_ops = [A1sgg, A1see, A2sgg, A2see]

results = mesolve(H, psi0, tspan, np.sqrt(Gamma) * A1sge, e_ops, options=Options(atol=1e-9,rtol=1e-9))

#plotting the results
fig, ax = plt.subplots(1,1)
ax.plot(tspan, results.expect[0], lw=2, label=r'$\sigma_{gg}^{A1}$')
ax.plot(tspan, results.expect[1], lw=2, label=r'$\sigma_{ee}^{A1}$')
ax.plot(tspan, results.expect[2], lw=2, label=r'$\sigma_{gg}^{A2}$')
ax.plot(tspan, results.expect[3], lw=2, label=r'$\sigma_{ee}^{A2}$')
ax.set_xlabel(r'Time $[\Omega^{-1}]$')
ax.legend()
plt.show()

Again changing the number of time steps in tspan to 500 gives the right result.

[nonhermitian]

I do not have much time to look into this until Monday. However, attached is a time-dependent three-level system, driven by a cavity , that works just fine. Maybe some clues can be found there.

#
# Single photon source based on a three level atom strongly coupled to a cavity
#
# We follow the treatment presented in Kuhn et al.,
# Appl. Phys. B 69, 373-377 (1999),
# http://www.mpq.mpg.de/qdynamics/publications/library/APB69p373_Kuhn.pdf,
# for more details see M. Hennrich's thesis,
# http://mediatum2.ub.tum.de/node?id=602970.
#
# We study the following lambda system,
#
#                |e>
#             --------
#             /     \
#      Omega /       \ g
#           /         \
#          /        -------
#      -------        |g>
#       |u>
#
# where |u> and |g> are the ground states and |e> is the exicted state.
# |u> and |e> are coupled by a classical laser field with Rabi frequency
# Omega, and |g> and |e> by a cavity field with 2g being the single-photon
# Rabi frequency.
#
from qutip import *
from pylab import *

# Define atomic states. Use ordering from paper
ustate = basis(3, 0)
excited = basis(3, 1)
ground = basis(3, 2)

# Set where to truncate Fock state for cavity
N = 2

# Create the atomic operators needed for the Hamiltonian
sigma_ge = tensor(qeye(N), ground * excited.dag())  # |g><e|
sigma_ue = tensor(qeye(N), ustate * excited.dag())  # |u><e|

# Create the photon operator
a = tensor(destroy(N), qeye(3))
ada = tensor(num(N), qeye(3))

# Define collapse operators
c_op_list = []
# Cavity decay rate
kappa = 1.5
c_op_list.append(sqrt(kappa) * a)

# Atomic decay rate
gamma = 6
# Use Rb branching ratio of 5/9 e->u, 4/9 e->g
c_op_list.append(sqrt(5 * gamma / 9) * sigma_ue)
c_op_list.append(sqrt(4 * gamma / 9) * sigma_ge)

# Define time vector
t = linspace(-15, 15, 100)
# Define pump strength as a function of time
wp = lambda t: 9 * exp(-(t / 5) ** 2)

# Set up the time varying Hamiltonian
g = 5
H0 = -g * (sigma_ge.dag() * a + a.dag() * sigma_ge)

H1 = (sigma_ue.dag() + sigma_ue)

H = [H0, [H1, '-9 * exp(-(t / 5) ** 2)']]

# Define initial state
psi0 = tensor(basis(N, 0), ustate)

# Define states onto which to project (same as in paper)
state_GG = tensor(basis(N, 1), ground)
sigma_GG = state_GG * state_GG.dag()
state_UU = tensor(basis(N, 0), ustate)
sigma_UU = state_UU * state_UU.dag()

output = mesolve(H, psi0, t, c_op_list,
                 [ada, sigma_UU, sigma_GG])

exp_ada, exp_uu, exp_gg = (output.expect[0], output.expect[1],
                           output.expect[2])

# Plot the results
figure()
subplot(211)
plot(t, wp(t), 'k')
ylabel('Control Field, $\Omega_\mathrm{p}$ [MHz]')
ax = twinx()
plot(t, kappa * exp_ada, 'b')
ylabel('Cavity emission rate, $1/\mu s$')
for tl in ax.get_yticklabels():
    tl.set_color('b')

subplot(212)
plot(t, exp_uu, 'k-', label='$P{\mathrm{uu}}$')
plot(t, exp_gg, 'k:', label='$P{\mathrm{gg}}$')
ylabel('Population')
xlabel('Time [$\mu s$]')
legend()
show()

[nlphysics]

Thanks. It's not that urgent. I'm just curious what the reason for that behavior might be. It's puzzling why does it apparently depend on the dimensions of dm as well as on the number of time steps.

[nonhermitian]

I quickly ran the problem, and it is still the slope is too large. Going from 100->50 in tanh resolves the issue. In order to look like a step pulse, you want the rate of change of the pulse to be much higher than the system time-scales. In your case, the system is on the time-scale of O(1), while, for tanh(100) the rate is ~60 at front and back, where as for tanh(50), it is ~30. In both cases, your system will see an effective step pulse.

My thinking is that if the argument of tanh is too large, then there is a big separation of time-scales in the problem. The ODE solver is trying to capture both the slow system dynamics, and the rapid step terms, and it is finding it difficult. This is a common issue with multiple, wide-ranging time scales.

On Oct 20, 2017, at 21:44, nlphysics notifications@github.com wrote:

Thanks. It's not that urgent. I'm just curious what the reason for that behavior might be. It's puzzling why does it apparently depend on the dimensions of dm as well as on the number of time steps.

— You are receiving this because you were assigned. Reply to this email directly, view it on GitHub https://github.com/qutip/qutip/issues/771#issuecomment-338361319, or mute the thread https://github.com/notifications/unsubscribe-auth/ABMPqbL-dF5jtatZLSmpxx_RYASR_2zEks5suWiugaJpZM4P__42.

[nonhermitian]

I am thinking that it could be that the solver is stepping over your pulse. Try setting the max_step parameter in options to something a fraction of the pulse width.

[nlphysics]

Your suspicion turned out to be right and setting the max_step parameter as you suggested solved the problem. Now it even works with a step pulse. Thank you for helping me and for all your effort!