Source code for qutip.floquet

# This file is part of QuTiP: Quantum Toolbox in Python.
#
#    Copyright (c) 2011 and later, Paul D. Nation and Robert J. Johansson.
#    All rights reserved.
#
#    Redistribution and use in source and binary forms, with or without
#    modification, are permitted provided that the following conditions are
#    met:
#
#    1. Redistributions of source code must retain the above copyright notice,
#       this list of conditions and the following disclaimer.
#
#    2. Redistributions in binary form must reproduce the above copyright
#       notice, this list of conditions and the following disclaimer in the
#       documentation and/or other materials provided with the distribution.
#
#    3. Neither the name of the QuTiP: Quantum Toolbox in Python nor the names
#       of its contributors may be used to endorse or promote products derived
#       from this software without specific prior written permission.
#
#    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
#    "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
#    LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
#    PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
#    HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
#    SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
#    LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
#    DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
#    THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
#    (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
#    OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
###############################################################################

__all__ = ['floquet_modes', 'floquet_modes_t', 'floquet_modes_table',
           'floquet_modes_t_lookup', 'floquet_states', 'floquet_states_t',
           'floquet_wavefunction', 'floquet_wavefunction_t',
           'floquet_state_decomposition', 'fsesolve',
           'floquet_master_equation_rates', 'floquet_collapse_operators',
           'floquet_master_equation_tensor',
           'floquet_master_equation_steadystate', 'floquet_basis_transform',
           'floquet_markov_mesolve', 'fmmesolve']

import numpy as np
import scipy.linalg as la
import scipy
from scipy import angle, pi, exp, sqrt
from types import FunctionType
from qutip.qobj import Qobj, isket
from qutip.superoperator import vec2mat_index, mat2vec, vec2mat
from qutip.mesolve import mesolve
from qutip.steadystate import steadystate
from qutip.states import ket2dm
from qutip.states import projection
from qutip.solver import Options
from qutip.propagator import propagator
from qutip.solver import Result
from qutip.cy.spmatfuncs import cy_ode_rhs
from qutip.expect import expect
from qutip.utilities import n_thermal


[docs]def floquet_modes(H, T, args=None, sort=False, U=None): """ Calculate the initial Floquet modes Phi_alpha(0) for a driven system with period T. Returns a list of :class:`qutip.qobj` instances representing the Floquet modes and a list of corresponding quasienergies, sorted by increasing quasienergy in the interval [-pi/T, pi/T]. The optional parameter `sort` decides if the output is to be sorted in increasing quasienergies or not. Parameters ---------- H : :class:`qutip.qobj` system Hamiltonian, time-dependent with period `T` args : dictionary dictionary with variables required to evaluate H T : float The period of the time-dependence of the hamiltonian. The default value 'None' indicates that the 'tlist' spans a single period of the driving. U : :class:`qutip.qobj` The propagator for the time-dependent Hamiltonian with period `T`. If U is `None` (default), it will be calculated from the Hamiltonian `H` using :func:`qutip.propagator.propagator`. Returns ------- output : list of kets, list of quasi energies Two lists: the Floquet modes as kets and the quasi energies. """ if U is None: # get the unitary propagator U = propagator(H, T, [], args) # find the eigenstates for the propagator evals, evecs = la.eig(U.full()) eargs = angle(evals) # make sure that the phase is in the interval [-pi, pi], so that # the quasi energy is in the interval [-pi/T, pi/T] where T is the # period of the driving. eargs += (eargs <= -2*pi) * (2*pi) + # (eargs > 0) * (-2*pi) eargs += (eargs <= -pi) * (2 * pi) + (eargs > pi) * (-2 * pi) e_quasi = -eargs / T # sort by the quasi energy if sort: order = np.argsort(-e_quasi) else: order = list(range(len(evals))) # prepare a list of kets for the floquet states new_dims = [U.dims[0], [1] * len(U.dims[0])] new_shape = [U.shape[0], 1] kets_order = [Qobj(np.matrix(evecs[:, o]).T, dims=new_dims, shape=new_shape) for o in order] return kets_order, e_quasi[order]
[docs]def floquet_modes_t(f_modes_0, f_energies, t, H, T, args=None): """ Calculate the Floquet modes at times tlist Phi_alpha(tlist) propagting the initial Floquet modes Phi_alpha(0) Parameters ---------- f_modes_0 : list of :class:`qutip.qobj` (kets) Floquet modes at :math:`t` f_energies : list Floquet energies. t : float The time at which to evaluate the floquet modes. H : :class:`qutip.qobj` system Hamiltonian, time-dependent with period `T` args : dictionary dictionary with variables required to evaluate H T : float The period of the time-dependence of the hamiltonian. Returns ------- output : list of kets The Floquet modes as kets at time :math:`t` """ # find t in [0,T] such that t_orig = t + n * T for integer n t = t - int(t / T) * T f_modes_t = [] # get the unitary propagator from 0 to t if t > 0.0: U = propagator(H, t, [], args) for n in np.arange(len(f_modes_0)): f_modes_t.append(U * f_modes_0[n] * exp(1j * f_energies[n] * t)) else: f_modes_t = f_modes_0 return f_modes_t
[docs]def floquet_modes_table(f_modes_0, f_energies, tlist, H, T, args=None): """ Pre-calculate the Floquet modes for a range of times spanning the floquet period. Can later be used as a table to look up the floquet modes for any time. Parameters ---------- f_modes_0 : list of :class:`qutip.qobj` (kets) Floquet modes at :math:`t` f_energies : list Floquet energies. tlist : array The list of times at which to evaluate the floquet modes. H : :class:`qutip.qobj` system Hamiltonian, time-dependent with period `T` T : float The period of the time-dependence of the hamiltonian. args : dictionary dictionary with variables required to evaluate H Returns ------- output : nested list A nested list of Floquet modes as kets for each time in `tlist` """ # truncate tlist to the driving period tlist_period = tlist[np.where(tlist <= T)] f_modes_table_t = [[] for t in tlist_period] opt = Options() opt.rhs_reuse = True for n, f_mode in enumerate(f_modes_0): output = mesolve(H, f_mode, tlist_period, [], [], args, opt) for t_idx, f_state_t in enumerate(output.states): f_modes_table_t[t_idx].append( f_state_t * exp(1j * f_energies[n] * tlist_period[t_idx])) return f_modes_table_t
[docs]def floquet_modes_t_lookup(f_modes_table_t, t, T): """ Lookup the floquet mode at time t in the pre-calculated table of floquet modes in the first period of the time-dependence. Parameters ---------- f_modes_table_t : nested list of :class:`qutip.qobj` (kets) A lookup-table of Floquet modes at times precalculated by :func:`qutip.floquet.floquet_modes_table`. t : float The time for which to evaluate the Floquet modes. T : float The period of the time-dependence of the hamiltonian. Returns ------- output : nested list A list of Floquet modes as kets for the time that most closely matching the time `t` in the supplied table of Floquet modes. """ # find t_wrap in [0,T] such that t = t_wrap + n * T for integer n t_wrap = t - int(t / T) * T # find the index in the table that corresponds to t_wrap (= tlist[t_idx]) t_idx = int(t_wrap / T * len(f_modes_table_t)) # XXX: might want to give a warning if the cast of t_idx to int discard # a significant fraction in t_idx, which would happen if the list of time # values isn't perfect matching the driving period # if debug: print "t = %f -> t_wrap = %f @ %d of %d" % (t, t_wrap, t_idx, # N) return f_modes_table_t[t_idx]
def floquet_states(f_modes_t, f_energies, t): """ Evaluate the floquet states at time t given the Floquet modes at that time. Parameters ---------- f_modes_t : list of :class:`qutip.qobj` (kets) A list of Floquet modes for time :math:`t`. f_energies : array The Floquet energies. t : float The time for which to evaluate the Floquet states. Returns ------- output : list A list of Floquet states for the time :math:`t`. """ return [(f_modes_t[i] * exp(-1j * f_energies[i] * t)) for i in np.arange(len(f_energies))]
[docs]def floquet_states_t(f_modes_0, f_energies, t, H, T, args=None): """ Evaluate the floquet states at time t given the initial Floquet modes. Parameters ---------- f_modes_t : list of :class:`qutip.qobj` (kets) A list of initial Floquet modes (for time :math:`t=0`). f_energies : array The Floquet energies. t : float The time for which to evaluate the Floquet states. H : :class:`qutip.qobj` System Hamiltonian, time-dependent with period `T`. T : float The period of the time-dependence of the hamiltonian. args : dictionary Dictionary with variables required to evaluate H. Returns ------- output : list A list of Floquet states for the time :math:`t`. """ f_modes_t = floquet_modes_t(f_modes_0, f_energies, t, H, T, args) return [(f_modes_t[i] * exp(-1j * f_energies[i] * t)) for i in np.arange(len(f_energies))]
def floquet_wavefunction(f_modes_t, f_energies, f_coeff, t): """ Evaluate the wavefunction for a time t using the Floquet state decompositon, given the Floquet modes at time `t`. Parameters ---------- f_modes_t : list of :class:`qutip.qobj` (kets) A list of initial Floquet modes (for time :math:`t=0`). f_energies : array The Floquet energies. f_coeff : array The coefficients for Floquet decomposition of the initial wavefunction. t : float The time for which to evaluate the Floquet states. Returns ------- output : :class:`qutip.qobj` The wavefunction for the time :math:`t`. """ return sum([f_modes_t[i] * exp(-1j * f_energies[i] * t) * f_coeff[i] for i in np.arange(len(f_energies))])
[docs]def floquet_wavefunction_t(f_modes_0, f_energies, f_coeff, t, H, T, args=None): """ Evaluate the wavefunction for a time t using the Floquet state decompositon, given the initial Floquet modes. Parameters ---------- f_modes_t : list of :class:`qutip.qobj` (kets) A list of initial Floquet modes (for time :math:`t=0`). f_energies : array The Floquet energies. f_coeff : array The coefficients for Floquet decomposition of the initial wavefunction. t : float The time for which to evaluate the Floquet states. H : :class:`qutip.qobj` System Hamiltonian, time-dependent with period `T`. T : float The period of the time-dependence of the hamiltonian. args : dictionary Dictionary with variables required to evaluate H. Returns ------- output : :class:`qutip.qobj` The wavefunction for the time :math:`t`. """ f_states_t = floquet_states_t(f_modes_0, f_energies, t, H, T, args) return sum([f_states_t[i] * f_coeff[i] for i in np.arange(len(f_energies))])
[docs]def floquet_state_decomposition(f_states, f_energies, psi): """ Decompose the wavefunction `psi` (typically an initial state) in terms of the Floquet states, :math:`\psi = \sum_\\alpha c_\\alpha \psi_\\alpha(0)`. Parameters ---------- f_states : list of :class:`qutip.qobj` (kets) A list of Floquet modes. f_energies : array The Floquet energies. psi : :class:`qutip.qobj` The wavefunction to decompose in the Floquet state basis. Returns ------- output : array The coefficients :math:`c_\\alpha` in the Floquet state decomposition. """ return [(f_states[i].dag() * psi).data[0, 0] for i in np.arange(len(f_energies))]
[docs]def fsesolve(H, psi0, tlist, e_ops=[], T=None, args={}, Tsteps=100): """ Solve the Schrodinger equation using the Floquet formalism. Parameters ---------- H : :class:`qutip.qobj.Qobj` System Hamiltonian, time-dependent with period `T`. psi0 : :class:`qutip.qobj` Initial state vector (ket). tlist : *list* / *array* list of times for :math:`t`. e_ops : list of :class:`qutip.qobj` / callback function list of operators for which to evaluate expectation values. If this list is empty, the state vectors for each time in `tlist` will be returned instead of expectation values. T : float The period of the time-dependence of the hamiltonian. args : dictionary Dictionary with variables required to evaluate H. Tsteps : integer The number of time steps in one driving period for which to precalculate the Floquet modes. `Tsteps` should be an even number. Returns ------- output : :class:`qutip.solver.Result` An instance of the class :class:`qutip.solver.Result`, which contains either an *array* of expectation values or an array of state vectors, for the times specified by `tlist`. """ if not T: # assume that tlist span exactly one period of the driving T = tlist[-1] # find the floquet modes for the time-dependent hamiltonian f_modes_0, f_energies = floquet_modes(H, T, args) # calculate the wavefunctions using the from the floquet modes f_modes_table_t = floquet_modes_table(f_modes_0, f_energies, np.linspace(0, T, Tsteps + 1), H, T, args) # setup Result for storing the results output = Result() output.times = tlist output.solver = "fsesolve" if isinstance(e_ops, FunctionType): output.num_expect = 0 expt_callback = True elif isinstance(e_ops, list): output.num_expect = len(e_ops) expt_callback = False if output.num_expect == 0: output.states = [] else: output.expect = [] for op in e_ops: if op.isherm: output.expect.append(np.zeros(len(tlist))) else: output.expect.append(np.zeros(len(tlist), dtype=complex)) else: raise TypeError("e_ops must be a list Qobj or a callback function") psi0_fb = psi0.transform(f_modes_0) for t_idx, t in enumerate(tlist): f_modes_t = floquet_modes_t_lookup(f_modes_table_t, t, T) f_states_t = floquet_states(f_modes_t, f_energies, t) psi_t = psi0_fb.transform(f_states_t, True) if expt_callback: # use callback method e_ops(t, psi_t) else: # calculate all the expectation values, or output psi if # no expectation value operators where defined if output.num_expect == 0: output.states.append(Qobj(psi_t)) else: for e_idx, e in enumerate(e_ops): output.expect[e_idx][t_idx] = expect(e, psi_t) return output
def floquet_master_equation_rates(f_modes_0, f_energies, c_op, H, T, args, J_cb, w_th, kmax=5, f_modes_table_t=None): """ Calculate the rates and matrix elements for the Floquet-Markov master equation. Parameters ---------- f_modes_0 : list of :class:`qutip.qobj` (kets) A list of initial Floquet modes. f_energies : array The Floquet energies. c_op : :class:`qutip.qobj` The collapse operators describing the dissipation. H : :class:`qutip.qobj` System Hamiltonian, time-dependent with period `T`. T : float The period of the time-dependence of the hamiltonian. args : dictionary Dictionary with variables required to evaluate H. J_cb : callback functions A callback function that computes the noise power spectrum, as a function of frequency, associated with the collapse operator `c_op`. w_th : float The temperature in units of frequency. k_max : int The truncation of the number of sidebands (default 5). f_modes_table_t : nested list of :class:`qutip.qobj` (kets) A lookup-table of Floquet modes at times precalculated by :func:`qutip.floquet.floquet_modes_table` (optional). Returns ------- output : list A list (Delta, X, Gamma, A) containing the matrices Delta, X, Gamma and A used in the construction of the Floquet-Markov master equation. """ N = len(f_energies) M = 2 * kmax + 1 omega = (2 * pi) / T Delta = np.zeros((N, N, M)) X = np.zeros((N, N, M), dtype=complex) Gamma = np.zeros((N, N, M)) A = np.zeros((N, N)) nT = 100 dT = T / nT tlist = np.arange(dT, T + dT / 2, dT) if f_modes_table_t is None: f_modes_table_t = floquet_modes_table(f_modes_0, f_energies, np.linspace(0, T, nT + 1), H, T, args) for t in tlist: # TODO: repeated invocations of floquet_modes_t is # inefficient... make a and b outer loops and use the mesolve # instead of the propagator. # f_modes_t = floquet_modes_t(f_modes_0, f_energies, t, H, T, args) f_modes_t = floquet_modes_t_lookup(f_modes_table_t, t, T) for a in range(N): for b in range(N): k_idx = 0 for k in range(-kmax, kmax + 1, 1): X[a, b, k_idx] += (dT / T) * exp(-1j * k * omega * t) * \ (f_modes_t[a].dag() * c_op * f_modes_t[b])[0, 0] k_idx += 1 Heaviside = lambda x: ((np.sign(x) + 1) / 2.0) for a in range(N): for b in range(N): k_idx = 0 for k in range(-kmax, kmax + 1, 1): Delta[a, b, k_idx] = f_energies[a] - f_energies[b] + k * omega Gamma[a, b, k_idx] = 2 * pi * Heaviside(Delta[a, b, k_idx]) * \ J_cb(Delta[a, b, k_idx]) * abs(X[a, b, k_idx]) ** 2 k_idx += 1 for a in range(N): for b in range(N): for k in range(-kmax, kmax + 1, 1): k1_idx = k + kmax k2_idx = -k + kmax A[a, b] += Gamma[a, b, k1_idx] + \ n_thermal(abs(Delta[a, b, k1_idx]), w_th) * \ (Gamma[a, b, k1_idx] + Gamma[b, a, k2_idx]) return Delta, X, Gamma, A def floquet_collapse_operators(A): """ Construct collapse operators corresponding to the Floquet-Markov master-equation rate matrix `A`. .. note:: Experimental. """ c_ops = [] N, M = np.shape(A) # # Here we really need a master equation on Bloch-Redfield form, or perhaps # we can use the Lindblad form master equation with some rotating frame # approximations? ... # for a in range(N): for b in range(N): if a != b and abs(A[a, b]) > 0.0: # only relaxation terms included... c_ops.append(sqrt(A[a, b]) * projection(N, a, b)) return c_ops def floquet_master_equation_tensor(Alist, f_energies): """ Construct a tensor that represents the master equation in the floquet basis (with constant Hamiltonian and collapse operators). Simplest RWA approximation [Grifoni et al, Phys.Rep. 304 229 (1998)] Parameters ---------- Alist : list A list of Floquet-Markov master equation rate matrices. f_energies : array The Floquet energies. Returns ------- output : array The Floquet-Markov master equation tensor `R`. """ if isinstance(Alist, list): # Alist can be a list of rate matrices corresponding # to different operators that couple to the environment N, M = np.shape(Alist[0]) else: # or a simple rate matrix, in which case we put it in a list Alist = [Alist] N, M = np.shape(Alist[0]) R = Qobj(scipy.sparse.csr_matrix((N * N, N * N)), [[N, N], [N, N]], [N * N, N * N]) R.data = R.data.tolil() for I in range(N * N): a, b = vec2mat_index(N, I) for J in range(N * N): c, d = vec2mat_index(N, J) R.data[I, J] = -1.0j * (f_energies[a] - f_energies[b]) * \ (a == c) * (b == d) for A in Alist: s1 = s2 = 0 for n in range(N): s1 += A[a, n] * (n == c) * (n == d) - A[n, a] * \ (a == c) * (a == d) s2 += (A[n, a] + A[n, b]) * (a == c) * (b == d) dR = (a == b) * s1 - 0.5 * (1 - (a == b)) * s2 if dR != 0.0: R.data[I, J] += dR R.data = R.data.tocsr() return R def floquet_master_equation_steadystate(H, A): """ Returns the steadystate density matrix (in the floquet basis!) for the Floquet-Markov master equation. """ c_ops = floquet_collapse_operators(A) rho_ss = steadystate(H, c_ops) return rho_ss def floquet_basis_transform(f_modes, f_energies, rho0): """ Make a basis transform that takes rho0 from the floquet basis to the computational basis. """ return rho0.transform(f_modes, True) # ----------------------------------------------------------------------------- # Floquet-Markov master equation # # def floquet_markov_mesolve(R, ekets, rho0, tlist, e_ops, f_modes_table=None, options=None, floquet_basis=True): """ Solve the dynamics for the system using the Floquet-Markov master equation. """ if options is None: opt = Options() else: opt = options if opt.tidy: R.tidyup() # # check initial state # if isket(rho0): # Got a wave function as initial state: convert to density matrix. rho0 = ket2dm(rho0) # # prepare output array # n_tsteps = len(tlist) dt = tlist[1] - tlist[0] output = Result() output.solver = "fmmesolve" output.times = tlist if isinstance(e_ops, FunctionType): n_expt_op = 0 expt_callback = True elif isinstance(e_ops, list): n_expt_op = len(e_ops) expt_callback = False if n_expt_op == 0: output.states = [] else: if not f_modes_table: raise TypeError("The Floquet mode table has to be provided " + "when requesting expectation values.") output.expect = [] output.num_expect = n_expt_op for op in e_ops: if op.isherm: output.expect.append(np.zeros(n_tsteps)) else: output.expect.append(np.zeros(n_tsteps, dtype=complex)) else: raise TypeError("Expectation parameter must be a list or a function") # # transform the initial density matrix to the eigenbasis: from # computational basis to the floquet basis # if ekets is not None: rho0 = rho0.transform(ekets) # # setup integrator # initial_vector = mat2vec(rho0.full()) r = scipy.integrate.ode(cy_ode_rhs) r.set_f_params(R.data.data, R.data.indices, R.data.indptr) r.set_integrator('zvode', method=opt.method, order=opt.order, atol=opt.atol, rtol=opt.rtol, max_step=opt.max_step) r.set_initial_value(initial_vector, tlist[0]) # # start evolution # rho = Qobj(rho0) t_idx = 0 for t in tlist: if not r.successful(): break rho.data = vec2mat(r.y) if expt_callback: # use callback method if floquet_basis: e_ops(t, Qobj(rho)) else: f_modes_table_t, T = f_modes_table f_modes_t = floquet_modes_t_lookup(f_modes_table_t, t, T) e_ops(t, Qobj(rho).transform(f_modes_t, True)) else: # calculate all the expectation values, or output rho if # no operators if n_expt_op == 0: if floquet_basis: output.states.append(Qobj(rho)) else: f_modes_table_t, T = f_modes_table f_modes_t = floquet_modes_t_lookup(f_modes_table_t, t, T) output.states.append(Qobj(rho).transform(f_modes_t, True)) else: f_modes_table_t, T = f_modes_table f_modes_t = floquet_modes_t_lookup(f_modes_table_t, t, T) for m in range(0, n_expt_op): output.expect[m][t_idx] = \ expect(e_ops[m], rho.transform(f_modes_t, False)) r.integrate(r.t + dt) t_idx += 1 return output # ----------------------------------------------------------------------------- # Solve the Floquet-Markov master equation # #
[docs]def fmmesolve(H, rho0, tlist, c_ops, e_ops=[], spectra_cb=[], T=None, args={}, options=Options(), floquet_basis=True, kmax=5): """ Solve the dynamics for the system using the Floquet-Markov master equation. .. note:: This solver currently does not support multiple collapse operators. Parameters ---------- H : :class:`qutip.qobj` system Hamiltonian. rho0 / psi0 : :class:`qutip.qobj` initial density matrix or state vector (ket). tlist : *list* / *array* list of times for :math:`t`. c_ops : list of :class:`qutip.qobj` list of collapse operators. e_ops : list of :class:`qutip.qobj` / callback function list of operators for which to evaluate expectation values. spectra_cb : list callback functions List of callback functions that compute the noise power spectrum as a function of frequency for the collapse operators in `c_ops`. T : float The period of the time-dependence of the hamiltonian. The default value 'None' indicates that the 'tlist' spans a single period of the driving. args : *dictionary* dictionary of parameters for time-dependent Hamiltonians and collapse operators. This dictionary should also contain an entry 'w_th', which is the temperature of the environment (if finite) in the energy/frequency units of the Hamiltonian. For example, if the Hamiltonian written in units of 2pi GHz, and the temperature is given in K, use the following conversion >>> temperature = 25e-3 # unit K >>> h = 6.626e-34 >>> kB = 1.38e-23 >>> args['w_th'] = temperature * (kB / h) * 2 * pi * 1e-9 options : :class:`qutip.solver` options for the ODE solver. k_max : int The truncation of the number of sidebands (default 5). Returns ------- output : :class:`qutip.solver` An instance of the class :class:`qutip.solver`, which contains either an *array* of expectation values for the times specified by `tlist`. """ if T is None: T = max(tlist) if len(spectra_cb) == 0: # add white noise callbacks if absent spectra_cb = [lambda w: 1.0] * len(c_ops) f_modes_0, f_energies = floquet_modes(H, T, args) f_modes_table_t = floquet_modes_table(f_modes_0, f_energies, np.linspace(0, T, 500 + 1), H, T, args) # get w_th from args if it exists if 'w_th' in args: w_th = args['w_th'] else: w_th = 0 # TODO: loop over input c_ops and spectra_cb, calculate one R for each set # calculate the rate-matrices for the floquet-markov master equation Delta, X, Gamma, Amat = floquet_master_equation_rates( f_modes_0, f_energies, c_ops[0], H, T, args, spectra_cb[0], w_th, kmax, f_modes_table_t) # the floquet-markov master equation tensor R = floquet_master_equation_tensor(Amat, f_energies) return floquet_markov_mesolve(R, f_modes_0, rho0, tlist, e_ops, f_modes_table=(f_modes_table_t, T), options=options, floquet_basis=floquet_basis)