FLiMESolve

doc/changes/2186.feature

@@ -0,0 +1 @@
+Added Floquet-Lindblad Master Equation for sinuisoidally driven time-dependent open system dynamics

doc/guide/dynamics/dynamics-floquet.rst

@@ -335,3 +335,54 @@ expectation values using:
     :context: reset
     :include-source: false
     :nofigs:
+
+The Floquet-Lindblad master equation in QuTiP
+---------------------------------------------
+The QuTiP function :func:`qutip.floquet.flimesolve` implements the 
+Floquet-Lindblad master equation. It calculates the dynamics of a system given 
+its initial state, a time-dependent Hamiltonian, and a list of operators 
+through which the system couples to its environment. 
+
+The following example extends the example studied above, and uses 
+:func:`qutip.floquet.flimesolve` to introduce dissipation into the calculation.
+
+
+.. plot:: guide/scripts/floquet_ex4.py
+   :width: 4.0in
+   :include-source:
+
+Importantly, the default solution method used by FLiMESolve will change, depending
+on the physical parameters of the system and the user inputs. A system
+with no time dependence can be solved in the full secular approximation, 
+whereas a system with complicated time-dependence will need a more relaxed 
+secular approximation. To adress this, FLiMESolve does have one additional input, 
+relative to MESolve. This property is called "time_sense," for time sensitivity,
+which allows for the secular approximation used in the Lindblad equation to 
+be relaxed. This is in contrast to FMMESolve, which uses the most restrictive 
+form of the secular approximation. The value of time sensitivity goes as
+
+.. math::
+    time sensitivity = (\omega_{1}-\omega_{2})/(S_{1}\S_{2})
+
+where :math:'\omega_{n},S_{n}' are the rotation frequency and rate of the 
+:math:'n^th' index of the Fourier-decomposed system operators in the Lindblad
+equation. Essentially, this value compares the rotation rate of a term in the
+master equation to how quickly that component affects the system. Large values
+mean that the component rotates quickly, relative to its contribution to the 
+system motion, so that overall the movement averages out on larger timescales.
+Smaller values have the converse meaning, such that overall the quotient can be
+understood to be a sort of "negligibility factor," with higher values being
+more negligible. The value of "time_sense" sets the value above which terms
+will be ignored. The default value of this input is zero, but if it can be
+to arbitrary limits to more accurately recover the behavior of MESolve.
+
+.. plot:: guide/scripts/floquet_ex5.py
+   :width: 4.0in
+   :include-source:
+
+Finally, for the sake of clarity, :func:`qutip.solver.floquet.fmmesolve`, 
+similar to :func:`qutip.solver.floquet.flimesolve`, always expects the 
+``e_ops`` to be specified in the laboratory basis:
+
+    output = flimesolve(H, psi0, tlist, [[sigmax() * gamma1**0.5 ]], e_ops=[num(2)], T=T, args=args)
+    p_ex = output.expect[0]

doc/guide/scripts/floquet_ex4.py

@@ -0,0 +1,57 @@
+import numpy as np
+from matplotlib import pyplot
+import qutip
+
+delta = 0.0 * 2 * np.pi
+eps0 = 1.0 * 2 * np.pi
+A = 0.25 * 2 * np.pi
+omega = 1.0 * 2 * np.pi
+T = 2 * np.pi / omega
+
+tlist = np.linspace(0.0, 20 * T, 301)
+psi0 = qutip.basis(2, 0)
+
+H0 = -delta / 2.0 * qutip.sigmax() - eps0 / 2.0 * qutip.sigmaz()
+H1 = A / 2.0 * qutip.sigmax()
+args = {"w": omega}
+H = [H0, [H1, lambda t, w: np.sin(w * t)]]
+gamma1 = 0.1
+
+
+# solve the floquet-lindblad master equation
+output = qutip.flimesolve(
+    H,
+    psi0,
+    tlist,
+    T,
+    c_ops=[np.sqrt(gamma1) * qutip.sigmax()],
+    args=args,
+    options={"store_floquet_states": True},
+)
+
+
+# calculate expectation values in the computational basis
+p_ex = np.zeros(tlist.shape, dtype=np.complex128)
+for idx, t in enumerate(tlist):
+    f_coeff_t = output.floquet_states[idx]
+    psi_t = output.floquet_basis.from_floquet_basis(f_coeff_t, t)
+    # Alternatively
+    psi_t = output.states[idx]
+    p_ex[idx] = qutip.expect(qutip.num(2), psi_t)
+
+# For reference: calculate the same thing with mesolve
+output = qutip.mesolve(
+    H, psi0, tlist, [np.sqrt(gamma1) * qutip.sigmax()], [qutip.num(2)], args
+)
+p_ex_ref = output.expect[0]
+
+
+# plot the results
+pyplot.plot(tlist, np.real(p_ex), "r--", tlist, 1 - np.real(p_ex), "b--")
+pyplot.plot(tlist, np.real(p_ex_ref), "r", tlist, 1 - np.real(p_ex_ref), "b")
+pyplot.xlabel("Time")
+pyplot.ylabel("Occupation probability")
+pyplot.legend(
+    ("Floquet $P_1$", "Floquet $P_0$", "Lindblad $P_1$", "Lindblad $P_0$")
+)
+pyplot.show()

doc/guide/scripts/floquet_ex5.py

@@ -0,0 +1,67 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Mon Apr 29 11:52:21 2024
+
+@author: frc00
+"""
+
+import numpy as np
+from matplotlib import pyplot
+import qutip
+
+delta = 0.0 * 2 * np.pi
+eps0 = 1.0 * 2 * np.pi
+A = 0.25 * 2 * np.pi
+omega = 1.0 * 2 * np.pi
+T = 2 * np.pi / omega
+
+tlist = np.linspace(0.0, 20 * T, 301)
+psi0 = qutip.basis(2, 0)
+
+H0 = -delta / 2.0 * qutip.sigmax() - eps0 / 2.0 * qutip.sigmaz()
+H1 = A / 2.0 * qutip.sigmax()
+args = {"w": omega}
+H = [H0, [H1, lambda t, w: np.sin(w * t)]]
+gamma1 = 0.1
+
+# setting the value of the time_sense argument arbitrarily high
+t_sensitivity = 1e10
+
+# solve the floquet-lindblad master equation
+output = qutip.flimesolve(
+    H,
+    psi0,
+    tlist,
+    T,
+    c_ops=[np.sqrt(gamma1) * qutip.sigmax()],
+    args=args,
+    options={"store_floquet_states": True},
+    time_sense=t_sensitivity,
+)
+
+
+# calculate expectation values in the computational basis
+p_ex = np.zeros(tlist.shape, dtype=np.complex128)
+for idx, t in enumerate(tlist):
+    f_coeff_t = output.floquet_states[idx]
+    psi_t = output.floquet_basis.from_floquet_basis(f_coeff_t, t)
+    # Alternatively
+    psi_t = output.states[idx]
+    p_ex[idx] = qutip.expect(qutip.num(2), psi_t)
+
+# For reference: calculate the same thing with mesolve
+output = qutip.mesolve(
+    H, psi0, tlist, [np.sqrt(gamma1) * qutip.sigmax()], [qutip.num(2)], args
+)
+p_ex_ref = output.expect[0]
+
+
+# plot the results
+pyplot.plot(tlist, np.real(p_ex), "r--", tlist, 1 - np.real(p_ex), "b--")
+pyplot.plot(tlist, np.real(p_ex_ref), "r", tlist, 1 - np.real(p_ex_ref), "b")
+pyplot.xlabel("Time")
+pyplot.ylabel("Occupation probability")
+pyplot.legend(
+    ("Floquet $P_1$", "Floquet $P_0$", "Lindblad $P_1$", "Lindblad $P_0$")
+)
+pyplot.show()

qutip/solver/__init__.py

@@ -11,6 +11,7 @@
 from .correlation import *
 from .spectrum import *
 from .floquet import *
+from .flimesolve import *
 from .floquet_bwcomp import *
 from .steadystate import *
 from .countstat import *