Propagator not working with time-dependent hamiltonian #2532
Description
Bug Description
Hi,
I've been trying to find the propagator for a time-dependent hamiltonian (with time-dependent collapse operators as well), and have been getting an error saying that the t_list lengths do not match in Qutip 5.0. I also tried to run it without using the collapse operators, but I get the same issue. The same code runs in Qutip 4.7.
My goal is to actually find two-time correlation functions of a time-dependent operator (following this paper:https://journals.aps.org/pra/abstract/10.1103/PhysRevA.102.023717). What I did do (and what has worked) is decompose the operator into time-independent parts and then use Qutip's correlation_2op_2t. This takes a really long time, since I need to do this for 8 different pairs of operators. I thought of instead finding the propagator and manually finding the correlation using it instead.
From the code below, you can see that mesolve works, but the propagator doesn't.
Code to Reproduce the Bug
gamma_1d_1 = Gamma
gamma_1d_2 = Gamma
## Input pulse
u = input_pulse
mod_u_2 = np.abs(u)**2
integral_u_2 = np.cumsum(mod_u_2) * t[1]
g_u = np.conjugate(u)/np.sqrt(1 - integral_u_2)
## Hamiltonian
emitter_dims = 3
mode_dims = 7
a = 1.1
psi_coherent = coherent(mode_dims, a)
## Mode dims is defined above
w0_1 = 2 * np.pi * 0
alpha_1 = 2 * np.pi * 0.3
w0_2 = 2 * np.pi * 0
alpha_2 = 2 * np.pi * 0.3
gamma_intrinsic = 2 * np.pi * 0.0
sig_1 = tensor(destroy(emitter_dims), qeye(emitter_dims), qeye(mode_dims))
sig_2 = tensor(qeye(emitter_dims), destroy(emitter_dims), qeye(mode_dims))
a_u = tensor(qeye(emitter_dims), qeye(emitter_dims), destroy(mode_dims))
H0 = (w0_1 * sig_1.dag() * sig_1 - alpha_1/2 * sig_1.dag() * sig_1.dag() * sig_1 * sig_1) + \
(w0_2 * sig_2.dag() * sig_2 - alpha_2/2 * sig_2.dag() * sig_2.dag() * sig_2 * sig_2) + \
0.5j * np.sqrt(gamma_1d_1 * gamma_1d_2) * (sig_1.dag() * sig_2 - sig_2.dag() * sig_1)
H_int1 = [0.5j * (a_u.dag() * (np.sqrt(gamma_1d_1) * sig_1 + np.sqrt(gamma_1d_2) * sig_2)), g_u]
H_int1_conj = [-0.5j * (a_u * (np.sqrt(gamma_1d_1) * sig_1.dag() + np.sqrt(gamma_1d_2) * sig_2.dag())), g_u.conjugate()]
gamma_phi_1 = 0.0 * gamma_1d_1
gamma_phi_2 = 0.0 * gamma_1d_2
Lo_us = [np.sqrt(gamma_1d_1)*sig_1, np.sqrt(gamma_1d_2)*sig_2, [a_u, g_u.conjugate()]]
# c_ops = [Lo_us, np.sqrt(gamma_phi_1) * sig_1.dag()*sig_1, np.sqrt(gamma_phi_2)*sig_2.dag()*sig_2]
c_ops = [Lo_us]
H = [H0, H_int1, H_int1_conj]
psi_0 = tensor(basis(emitter_dims, 0), basis(emitter_dims, 0), psi_coherent)
results = mesolve(H, psi_0, t, c_ops = c_ops)
ρ0 = psi_0 * psi_0.dag()
U = propagator(H, t, c_ops=c_ops)Code Output
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
Cell In[10], line 2
1 ρ0 = psi_0 * psi_0.dag()
----> 2 U = propagator(H, t, c_ops=c_ops)
File /opt/anaconda3/envs/qutip5/lib/python3.11/site-packages/qutip/solver/propagator.py:66, in propagator(H, t, c_ops, args, options, **kwargs)
63 list_output = True
65 if not isinstance(H, (Qobj, QobjEvo)):
---> 66 H = QobjEvo(H, args=args, **kwargs)
68 if c_ops:
69 H = liouvillian(H, c_ops)
File /opt/anaconda3/envs/qutip5/lib/python3.11/site-packages/qutip/core/cy/qobjevo.pyx:223, in qutip.core.cy.qobjevo.QobjEvo.__init__()
File /opt/anaconda3/envs/qutip5/lib/python3.11/site-packages/qutip/core/cy/qobjevo.pyx:272, in qutip.core.cy.qobjevo.QobjEvo._read_element()
File /opt/anaconda3/envs/qutip5/lib/python3.11/site-packages/qutip/core/coefficient.py:175, in coefficient(base, tlist, args, args_ctypes, order, compile_opt, function_style, boundary_conditions, **kwargs)
173 for type_ in coefficient_builders:
174 if isinstance(base, type_):
--> 175 return coefficient_builders[type_](base, **kwargs)
177 if callable(base):
178 op = FunctionCoefficient(base, args.copy(), style=function_style)
File /opt/anaconda3/envs/qutip5/lib/python3.11/site-packages/qutip/core/cy/coefficient.pyx:418, in qutip.core.cy.coefficient.InterCoefficient.__init__()
ValueError: tlist must be the same len as the array to interpolateExpected Behaviour
Wanted it to return the propagator
Your Environment
QuTiP Version: 5.0.4
Numpy Version: 1.26.4
Scipy Version: 1.12.0
Cython Version: None
Matplotlib Version: 3.9.2
Python Version: 3.11.9
Number of CPUs: 8
BLAS Info: Generic
INTEL MKL Ext: False
Platform Info: Darwin (x86_64)Additional Context
Even in Qutip 4.7, though the code runs, I'm seeing quite different outputs from mesolve and using the propagator (I might be using the propagator incorrectly though, so some help is appreciated):
plt.plot(t, expect(a_u.dag() * a_u, results.states), label="mesolve")
for i in range(len(t)):
plt.plot(t[i], expect(a_u.dag() * a_u, U[i]*ρ0), "x", color="orange")
plt.ylabel("Population")
plt.xlabel("Time")
plt.legend()
plt.show()Output: