I am trying to implement some codes to describe the coherent control of a system. The problem is that I want to use an hamiltonian of the form: {H}_ij = u(t) * h_ij * exp(omega_ij * t) for the i^th row and j^th column with h_ij = conj(h_ji) and h_ii = 0. Because the time dependence of each terms is different, this hamiltonian cannot be written using the techniques to write a time-dependent hamiltonian with the first three options detailed in the user manual (http://qutip.org/downloads/4.4.0/qutip-doc-4.4.pdf#subsection.3.5.6) more precisely at page 61.
It can however be realized with the fourth (outdated) option:
Finally option #4, expressing the Hamiltonian as a Python function, is the original method for time dependence in QuTiP 1.x. However, this method is somewhat less efficient then the previously mentioned methods. However, in contrast to the other options this method can be used in implementing time-dependent Hamiltonians that cannot be expressed as a function of constant operators with time-dependent coefficients.
I think this is exactly what I want but I don't understand how I can implement the Hamiltonian as a python function. Can someone show me an example of such "Hamiltonian as a Python function"?
EDIT: This is what I've tried so far:
import numpy as np
import qutip
psi_0 = qutip.ket2dm(qutip.basis(2, 0))
psi_1 = qutip.ket2dm(qutip.basis(2, 1))
time = np.linspace(0, 10)
H_zeros = qutip.Qobj([[0, 0], [0, 0]])
def Hamiltonian_1(t, args):
return [[0, np.exp(1j * t)], [np.exp(- 1j * t), 0]]
def Hamiltonian_2(t):
return [[0, np.exp(1j * t)], [np.exp(- 1j * t), 0]]
def Hamiltonian_3(t, args):
return qutip.Qobj([[0, np.exp(1j * t)], [np.exp(- 1j * t), 0]])
def Hamiltonian_4(t):
return qutip.Qobj([[0, np.exp(1j * t)], [np.exp(- 1j * t), 0]])
# result = qutip.mesolve(Hamiltonian_1, psi_0, time, [], [psi_0, psi_1])
# result = qutip.mesolve(Hamiltonian_2, psi_0, time, [], [psi_0, psi_1])
# result = qutip.mesolve(Hamiltonian_3, psi_0, time, [], [psi_0, psi_1])
# result = qutip.mesolve(Hamiltonian_4, psi_0, time, [], [psi_0, psi_1])
# result = qutip.mesolve([Hamiltonian_1], psi_0, time, [], [psi_0, psi_1])
# result = qutip.mesolve([Hamiltonian_2], psi_0, time, [], [psi_0, psi_1])
# result = qutip.mesolve([Hamiltonian_3], psi_0, time, [], [psi_0, psi_1])
# result = qutip.mesolve([Hamiltonian_4], psi_0, time, [], [psi_0, psi_1])
# result = qutip.mesolve([H_zeros, Hamiltonian_1], psi_0, time, [], [psi_0, psi_1])
# result = qutip.mesolve([H_zeros, Hamiltonian_2], psi_0, time, [], [psi_0, psi_1])
# result = qutip.mesolve([H_zeros, Hamiltonian_3], psi_0, time, [], [psi_0, psi_1])
# result = qutip.mesolve([H_zeros, Hamiltonian_4], psi_0, time, [], [psi_0, psi_1])
# result = qutip.mesolve([H_zeros, [H_zeros, Hamiltonian_1]], psi_0, time, [], [psi_0, psi_1])
# result = qutip.mesolve([H_zeros, [H_zeros, Hamiltonian_2]], psi_0, time, [], [psi_0, psi_1])
# result = qutip.mesolve([H_zeros, [H_zeros, Hamiltonian_3]], psi_0, time, [], [psi_0, psi_1])
# result = qutip.mesolve([H_zeros, [H_zeros, Hamiltonian_4]], psi_0, time, [], [psi_0, psi_1])
I've tried in this simple case with different definitions of the hamiltonian function (Hamiltonian_1 to 4) and I've tried multiple combination of mesolve functions (commented in this code for the moment)... None of them is working.
Thank you, Paul