diff --git a/doc/QuTiP_tree_plot/d3_data/qutip.json b/doc/QuTiP_tree_plot/d3_data/qutip.json
index f93b96f9bd..0f46d672ee 100644
--- a/doc/QuTiP_tree_plot/d3_data/qutip.json
+++ b/doc/QuTiP_tree_plot/d3_data/qutip.json
@@ -1 +1 @@
-{"name": "QuTiP", "color": "black", "index": -1, "children": [{"name": "core", "color": "black", "index": -1, "children": [{"name": "cy", "color": "black", "index": -1, "children": [{"name": "coefficient", "color": "black", "index": -1, "children": [{"name": "Coefficient", "color": "black", "index": -1}, {"name": "ConjCoefficient", "color": "black", "index": -1}, {"name": "FunctionCoefficient", "color": "black", "index": -1}, {"name": "InterCoefficient", "color": "black", "index": -1}, {"name": "MulCoefficient", "color": "black", "index": -1}, {"name": "NormCoefficient", "color": "black", "index": -1}, {"name": "ShiftCoefficient", "color": "black", "index": -1}, {"name": "StrFunctionCoefficient", "color": "black", "index": -1}, {"name": "SumCoefficient", "color": "black", "index": -1}, {"name": "coefficient_function_parameters", "color": "black", "index": -1}, {"name": "proj", "color": "black", "index": -1}]}, {"name": "_element", "color": "black", "index": -1}, {"name": "qobjevo", "color": "black", "index": -1, "children": [{"name": "QobjEvo", "color": "black", "index": -1}]}]}, {"name": "data", "color": "black", "index": -1, "children": [{"name": "csr", "color": "black", "index": -1, "children": [{"name": "CSR", "color": "black", "index": -1}, {"name": "Sorter", "color": "black", "index": -1}, {"name": "copy_structure", "color": "black", "index": -1}, {"name": "diags", "color": "black", "index": -1}, {"name": "empty", "color": "black", "index": -1}, {"name": "empty_like", "color": "black", "index": -1}, {"name": "fast_from_scipy", "color": "black", "index": -1}, {"name": "from_dense", "color": "black", "index": -1}, {"name": "identity", "color": "black", "index": -1}, {"name": "nnz", "color": "black", "index": -1}, {"name": "sorted", "color": "black", "index": -1}, {"name": "zeros", "color": "black", "index": -1}]}, {"name": "base", "color": "black", "index": -1, "children": [{"name": "Data", "color": "black", "index": -1}, {"name": "EfficiencyWarning", "color": "black", "index": -1}]}, {"name": "dense", "color": "black", "index": -1, "children": [{"name": "Dense", "color": "black", "index": -1}, {"name": "OrderEfficiencyWarning", "color": "black", "index": -1}, {"name": "diags", "color": "black", "index": -1}, {"name": "empty", "color": "black", "index": -1}, {"name": "empty_like", "color": "black", "index": -1}, {"name": "fast_from_numpy", "color": "black", "index": -1}, {"name": "from_csr", "color": "black", "index": -1}, {"name": "identity", "color": "black", "index": -1}, {"name": "zeros", "color": "black", "index": -1}]}, {"name": "dispatch", "color": "black", "index": -1, "children": [{"name": "Dispatcher", "color": "black", "index": -1}]}, {"name": "add", "color": "black", "index": -1, "children": [{"name": "add_csr", "color": "black", "index": -1}, {"name": "add_dense", "color": "black", "index": -1}, {"name": "iadd_dense", "color": "black", "index": -1}, {"name": "sub_csr", "color": "black", "index": -1}, {"name": "sub_dense", "color": "black", "index": -1}]}, {"name": "adjoint", "color": "black", "index": -1, "children": [{"name": "adjoint_csr", "color": "black", "index": -1}, {"name": "adjoint_dense", "color": "black", "index": -1}, {"name": "conj_csr", "color": "black", "index": -1}, {"name": "conj_dense", "color": "black", "index": -1}, {"name": "transpose_csr", "color": "black", "index": -1}, {"name": "transpose_dense", "color": "black", "index": -1}]}, {"name": "reshape", "color": "black", "index": -1, "children": [{"name": "column_stack_csr", "color": "black", "index": -1}, {"name": "column_stack_dense", "color": "black", "index": -1}, {"name": "column_unstack_csr", "color": "black", "index": -1}, {"name": "column_unstack_dense", "color": "black", "index": -1}, {"name": "reshape_csr", "color": "black", "index": -1}, {"name": "reshape_dense", "color": "black", "index": -1}, {"name": "split_columns_csr", "color": "black", "index": -1}, {"name": "split_columns_dense", "color": "black", "index": -1}]}, {"name": "constant", "color": "black", "index": -1}, {"name": "convert", "color": "black", "index": -1}, {"name": "eigen", "color": "black", "index": -1, "children": [{"name": "eigs_csr", "color": "black", "index": -1}, {"name": "eigs_dense", "color": "black", "index": -1}]}, {"name": "expect", "color": "black", "index": -1, "children": [{"name": "expect_csr", "color": "black", "index": -1}, {"name": "expect_csr_dense", "color": "black", "index": -1}, {"name": "expect_dense", "color": "black", "index": -1}, {"name": "expect_super_csr", "color": "black", "index": -1}, {"name": "expect_super_csr_dense", "color": "black", "index": -1}, {"name": "expect_super_dense", "color": "black", "index": -1}]}, {"name": "expm", "color": "black", "index": -1, "children": [{"name": "expm_csr", "color": "black", "index": -1}, {"name": "expm_csr_dense", "color": "black", "index": -1}]}, {"name": "properties", "color": "black", "index": -1, "children": [{"name": "isdiag_csr", "color": "black", "index": -1}, {"name": "isherm_csr", "color": "black", "index": -1}, {"name": "iszero_csr", "color": "black", "index": -1}, {"name": "iszero_dense", "color": "black", "index": -1}]}, {"name": "mul", "color": "black", "index": -1, "children": [{"name": "imul_csr", "color": "black", "index": -1}, {"name": "imul_data", "color": "black", "index": -1}, {"name": "imul_dense", "color": "black", "index": -1}, {"name": "mul_csr", "color": "black", "index": -1}, {"name": "mul_dense", "color": "black", "index": -1}, {"name": "neg_csr", "color": "black", "index": -1}, {"name": "neg_dense", "color": "black", "index": -1}]}, {"name": "inner", "color": "black", "index": -1, "children": [{"name": "inner_csr", "color": "black", "index": -1}, {"name": "inner_op_csr", "color": "black", "index": -1}]}, {"name": "linalg", "color": "black", "index": -1, "children": [{"name": "inv_csr", "color": "black", "index": -1}, {"name": "inv_dense", "color": "black", "index": -1}]}, {"name": "kron", "color": "black", "index": -1, "children": [{"name": "kron_csr", "color": "black", "index": -1}, {"name": "kron_dense", "color": "black", "index": -1}]}, {"name": "make", "color": "black", "index": -1, "children": [{"name": "one_element_csr", "color": "black", "index": -1}, {"name": "one_element_dense", "color": "black", "index": -1}]}, {"name": "matmul", "color": "black", "index": -1, "children": [{"name": "matmul_csr", "color": "black", "index": -1}, {"name": "matmul_csr_dense_dense", "color": "black", "index": -1}, {"name": "matmul_dense", "color": "black", "index": -1}, {"name": "multiply_csr", "color": "black", "index": -1}, {"name": "multiply_dense", "color": "black", "index": -1}]}, {"name": "norm", "color": "black", "index": -1, "children": [{"name": "frobenius_csr", "color": "black", "index": -1}, {"name": "frobenius_data", "color": "black", "index": -1}, {"name": "frobenius_dense", "color": "black", "index": -1}, {"name": "l2_csr", "color": "black", "index": -1}, {"name": "l2_dense", "color": "black", "index": -1}, {"name": "max_csr", "color": "black", "index": -1}, {"name": "max_dense", "color": "black", "index": -1}, {"name": "one_csr", "color": "black", "index": -1}, {"name": "one_dense", "color": "black", "index": -1}, {"name": "trace_csr", "color": "black", "index": -1}, {"name": "trace_dense", "color": "black", "index": -1}]}, {"name": "permute", "color": "black", "index": -1, "children": [{"name": "dimensions_csr", "color": "black", "index": -1}, {"name": "indices_csr", "color": "black", "index": -1}]}, {"name": "pow", "color": "black", "index": -1, "children": [{"name": "pow_csr", "color": "black", "index": -1}]}, {"name": "project", "color": "black", "index": -1, "children": [{"name": "project_csr", "color": "black", "index": -1}, {"name": "project_dense", "color": "black", "index": -1}]}, {"name": "ptrace", "color": "black", "index": -1, "children": [{"name": "ptrace_csr", "color": "black", "index": -1}, {"name": "ptrace_csr_dense", "color": "black", "index": -1}, {"name": "ptrace_dense", "color": "black", "index": -1}]}, {"name": "tidyup", "color": "black", "index": -1, "children": [{"name": "tidyup_csr", "color": "black", "index": -1}, {"name": "tidyup_dense", "color": "black", "index": -1}]}, {"name": "trace", "color": "black", "index": -1, "children": [{"name": "trace_csr", "color": "black", "index": -1}, {"name": "trace_dense", "color": "black", "index": -1}]}]}, {"name": "coefficient", "color": "black", "index": -1, "children": [{"name": "CompilationOptions", "color": "black", "index": -1}, {"name": "StringParsingWarning", "color": "black", "index": -1}, {"name": "clean_compiled_coefficient", "color": "black", "index": -1}, {"name": "coeff_from_str", "color": "black", "index": -1}, {"name": "coefficient", "color": "black", "index": -1}, {"name": "compileType", "color": "black", "index": -1}, {"name": "compile_code", "color": "black", "index": -1}, {"name": "conj", "color": "black", "index": -1}, {"name": "extract_constant", "color": "black", "index": -1}, {"name": "extract_cte_pattern", "color": "black", "index": -1}, {"name": "find_type_from_str", "color": "black", "index": -1}, {"name": "fix_type", "color": "black", "index": -1}, {"name": "fromstr", "color": "black", "index": -1}, {"name": "make_cy_code", "color": "black", "index": -1}, {"name": "norm", "color": "black", "index": -1}, {"name": "parse", "color": "black", "index": -1}, {"name": "proj", "color": "black", "index": -1}, {"name": "shift", "color": "black", "index": -1}, {"name": "space_parts", "color": "black", "index": -1}, {"name": "test_parsed", "color": "black", "index": -1}, {"name": "try_import", "color": "black", "index": -1}, {"name": "try_parse", "color": "black", "index": -1}, {"name": "use_hinted_type", "color": "black", "index": -1}]}, {"name": "options", "color": "black", "index": -1, "children": [{"name": "CoreOptions", "color": "black", "index": -1}, {"name": "QutipOptions", "color": "black", "index": -1}]}, {"name": "qobj", "color": "black", "index": -1, "children": [{"name": "Qobj", "color": "black", "index": -1}, {"name": "isbra", "color": "black", "index": -1}, {"name": "isherm", "color": "black", "index": -1}, {"name": "isket", "color": "black", "index": -1}, {"name": "isoper", "color": "black", "index": -1}, {"name": "isoperbra", "color": "black", "index": -1}, {"name": "isoperket", "color": "black", "index": -1}, {"name": "issuper", "color": "black", "index": -1}, {"name": "ptrace", "color": "black", "index": -1}]}, {"name": "dimensions", "color": "black", "index": -1, "children": [{"name": "collapse_dims_oper", "color": "black", "index": -1}, {"name": "collapse_dims_super", "color": "black", "index": -1}, {"name": "deep_map", "color": "black", "index": -1}, {"name": "deep_remove", "color": "black", "index": -1}, {"name": "dims_idxs_to_tensor_idxs", "color": "black", "index": -1}, {"name": "dims_to_tensor_perm", "color": "black", "index": -1}, {"name": "dims_to_tensor_shape", "color": "black", "index": -1}, {"name": "enumerate_flat", "color": "black", "index": -1}, {"name": "flatten", "color": "black", "index": -1}, {"name": "is_scalar", "color": "black", "index": -1}, {"name": "is_vector", "color": "black", "index": -1}, {"name": "is_vectorized_oper", "color": "black", "index": -1}, {"name": "type_from_dims", "color": "black", "index": -1}, {"name": "unflatten", "color": "black", "index": -1}]}, {"name": "metrics", "color": "black", "index": -1, "children": [{"name": "average_gate_fidelity", "color": "black", "index": -1}, {"name": "bures_angle", "color": "black", "index": -1}, {"name": "bures_dist", "color": "black", "index": -1}, {"name": "dnorm", "color": "black", "index": -1}, {"name": "fidelity", "color": "black", "index": -1}, {"name": "hellinger_dist", "color": "black", "index": -1}, {"name": "hilbert_dist", "color": "black", "index": -1}, {"name": "process_fidelity", "color": "black", "index": -1}, {"name": "tracedist", "color": "black", "index": -1}, {"name": "unitarity", "color": "black", "index": -1}]}, {"name": "superop_reps", "color": "black", "index": -1, "children": [{"name": "isqubitdims", "color": "black", "index": -1}, {"name": "kraus_to_choi", "color": "black", "index": -1}, {"name": "kraus_to_super", "color": "black", "index": -1}, {"name": "to_chi", "color": "black", "index": -1}, {"name": "to_choi", "color": "black", "index": -1}, {"name": "to_kraus", "color": "black", "index": -1}, {"name": "to_stinespring", "color": "black", "index": -1}, {"name": "to_super", "color": "black", "index": -1}]}, {"name": "states", "color": "black", "index": -1, "children": [{"name": "basis", "color": "black", "index": -1}, {"name": "bell_state", "color": "black", "index": -1}, {"name": "bra", "color": "black", "index": -1}, {"name": "coherent", "color": "black", "index": -1}, {"name": "coherent_dm", "color": "black", "index": -1}, {"name": "enr_fock", "color": "black", "index": -1}, {"name": "enr_state_dictionaries", "color": "black", "index": -1}, {"name": "enr_thermal_dm", "color": "black", "index": -1}, {"name": "fock", "color": "black", "index": -1}, {"name": "fock_dm", "color": "black", "index": -1}, {"name": "ghz_state", "color": "black", "index": -1}, {"name": "ket", "color": "black", "index": -1}, {"name": "ket2dm", "color": "black", "index": -1}, {"name": "maximally_mixed_dm", "color": "black", "index": -1}, {"name": "phase_basis", "color": "black", "index": -1}, {"name": "projection", "color": "black", "index": -1}, {"name": "qstate", "color": "black", "index": -1}, {"name": "qutrit_basis", "color": "black", "index": -1}, {"name": "singlet_state", "color": "black", "index": -1}, {"name": "spin_coherent", "color": "black", "index": -1}, {"name": "spin_state", "color": "black", "index": -1}, {"name": "state_index_number", "color": "black", "index": -1}, {"name": "state_number_enumerate", "color": "black", "index": -1}, {"name": "state_number_index", "color": "black", "index": -1}, {"name": "state_number_qobj", "color": "black", "index": -1}, {"name": "thermal_dm", "color": "black", "index": -1}, {"name": "triplet_states", "color": "black", "index": -1}, {"name": "w_state", "color": "black", "index": -1}, {"name": "zero_ket", "color": "black", "index": -1}]}, {"name": "operators", "color": "black", "index": -1, "children": [{"name": "charge", "color": "black", "index": -1}, {"name": "commutator", "color": "black", "index": -1}, {"name": "create", "color": "black", "index": -1}, {"name": "destroy", "color": "black", "index": -1}, {"name": "displace", "color": "black", "index": -1}, {"name": "enr_destroy", "color": "black", "index": -1}, {"name": "enr_identity", "color": "black", "index": -1}, {"name": "qeye", "color": "black", "index": -1}, {"name": "jmat", "color": "black", "index": -1}, {"name": "momentum", "color": "black", "index": -1}, {"name": "num", "color": "black", "index": -1}, {"name": "phase", "color": "black", "index": -1}, {"name": "position", "color": "black", "index": -1}, {"name": "qdiags", "color": "black", "index": -1}, {"name": "qft", "color": "black", "index": -1}, {"name": "qutrit_ops", "color": "black", "index": -1}, {"name": "qzero", "color": "black", "index": -1}, {"name": "sigmam", "color": "black", "index": -1}, {"name": "sigmap", "color": "black", "index": -1}, {"name": "sigmax", "color": "black", "index": -1}, {"name": "sigmay", "color": "black", "index": -1}, {"name": "sigmaz", "color": "black", "index": -1}, {"name": "spin_J_set", "color": "black", "index": -1}, {"name": "spin_Jm", "color": "black", "index": -1}, {"name": "spin_Jp", "color": "black", "index": -1}, {"name": "spin_Jx", "color": "black", "index": -1}, {"name": "spin_Jy", "color": "black", "index": -1}, {"name": "spin_Jz", "color": "black", "index": -1}, {"name": "squeeze", "color": "black", "index": -1}, {"name": "squeezing", "color": "black", "index": -1}, {"name": "tunneling", "color": "black", "index": -1}]}, {"name": "tensor", "color": "black", "index": -1, "children": [{"name": "composite", "color": "black", "index": -1}, {"name": "expand_operator", "color": "black", "index": -1}, {"name": "super_tensor", "color": "black", "index": -1}, {"name": "tensor", "color": "black", "index": -1}, {"name": "tensor_contract", "color": "black", "index": -1}, {"name": "tensor_swap", "color": "black", "index": -1}]}, {"name": "superoperator", "color": "black", "index": -1, "children": [{"name": "lindblad_dissipator", "color": "black", "index": -1}, {"name": "liouvillian", "color": "black", "index": -1}, {"name": "operator_to_vector", "color": "black", "index": -1}, {"name": "reshuffle", "color": "black", "index": -1}, {"name": "spost", "color": "black", "index": -1}, {"name": "spre", "color": "black", "index": -1}, {"name": "sprepost", "color": "black", "index": -1}, {"name": "stack_columns", "color": "black", "index": -1}, {"name": "stacked_index", "color": "black", "index": -1}, {"name": "unstack_columns", "color": "black", "index": -1}, {"name": "unstacked_index", "color": "black", "index": -1}, {"name": "vector_to_operator", "color": "black", "index": -1}]}, {"name": "semidefinite", "color": "black", "index": -1, "children": [{"name": "Complex", "color": "black", "index": -1}, {"name": "bmat", "color": "black", "index": -1}, {"name": "complex_var", "color": "black", "index": -1}, {"name": "conj", "color": "black", "index": -1}, {"name": "dag", "color": "black", "index": -1}, {"name": "dens", "color": "black", "index": -1}, {"name": "dnorm_problem", "color": "black", "index": -1}, {"name": "dnorm_sparse_problem", "color": "black", "index": -1}, {"name": "herm", "color": "black", "index": -1}, {"name": "initialize_constraints_on_dnorm_problem", "color": "black", "index": -1}, {"name": "kron", "color": "black", "index": -1}, {"name": "memoize", "color": "black", "index": -1}, {"name": "pos", "color": "black", "index": -1}, {"name": "pos_noherm", "color": "black", "index": -1}, {"name": "qudit_swap", "color": "black", "index": -1}]}, {"name": "_brtensor", "color": "black", "index": -1}, {"name": "_brtools", "color": "black", "index": -1, "children": [{"name": "SpectraCoefficient", "color": "black", "index": -1}, {"name": "matmul_var_data", "color": "black", "index": -1}]}, {"name": "blochredfield", "color": "black", "index": -1, "children": [{"name": "bloch_redfield_tensor", "color": "black", "index": -1}, {"name": "brterm", "color": "black", "index": -1}]}, {"name": "expect", "color": "black", "index": -1, "children": [{"name": "expect", "color": "black", "index": -1}, {"name": "variance", "color": "black", "index": -1}]}, {"name": "gates", "color": "black", "index": -1, "children": [{"name": "berkeley", "color": "black", "index": -1}, {"name": "cnot", "color": "black", "index": -1}, {"name": "cphase", "color": "black", "index": -1}, {"name": "cs_gate", "color": "black", "index": -1}, {"name": "csign", "color": "black", "index": -1}, {"name": "ct_gate", "color": "black", "index": -1}, {"name": "cy_gate", "color": "black", "index": -1}, {"name": "cz_gate", "color": "black", "index": -1}, {"name": "fredkin", "color": "black", "index": -1}, {"name": "globalphase", "color": "black", "index": -1}, {"name": "hadamard_transform", "color": "black", "index": -1}, {"name": "iswap", "color": "black", "index": -1}, {"name": "molmer_sorensen", "color": "black", "index": -1}, {"name": "phasegate", "color": "black", "index": -1}, {"name": "qrot", "color": "black", "index": -1}, {"name": "qubit_clifford_group", "color": "black", "index": -1}, {"name": "rx", "color": "black", "index": -1}, {"name": "ry", "color": "black", "index": -1}, {"name": "rz", "color": "black", "index": -1}, {"name": "s_gate", "color": "black", "index": -1}, {"name": "snot", "color": "black", "index": -1}, {"name": "sqrtiswap", "color": "black", "index": -1}, {"name": "sqrtnot", "color": "black", "index": -1}, {"name": "sqrtswap", "color": "black", "index": -1}, {"name": "swap", "color": "black", "index": -1}, {"name": "swapalpha", "color": "black", "index": -1}, {"name": "t_gate", "color": "black", "index": -1}, {"name": "toffoli", "color": "black", "index": -1}]}, {"name": "subsystem_apply", "color": "black", "index": -1, "children": [{"name": "subsystem_apply", "color": "black", "index": -1}]}]}, {"name": "optionsclass", "color": "black", "index": -1, "children": [{"name": "AllOptions", "color": "black", "index": -1}, {"name": "optionsclass", "color": "black", "index": -1}]}, {"name": "solve", "color": "black", "index": -1, "children": [{"name": "piqs", "color": "black", "index": -1, "children": [{"name": "Dicke", "color": "black", "index": -1}, {"name": "Pim", "color": "black", "index": -1}, {"name": "am", "color": "black", "index": -1}, {"name": "ap", "color": "black", "index": -1}, {"name": "block_matrix", "color": "black", "index": -1}, {"name": "collapse_uncoupled", "color": "black", "index": -1}, {"name": "css", "color": "black", "index": -1}, {"name": "dicke", "color": "black", "index": -1}, {"name": "dicke_basis", "color": "black", "index": -1}, {"name": "dicke_blocks", "color": "black", "index": -1}, {"name": "dicke_blocks_full", "color": "black", "index": -1}, {"name": "dicke_function_trace", "color": "black", "index": -1}, {"name": "energy_degeneracy", "color": "black", "index": -1}, {"name": "entropy_vn_dicke", "color": "black", "index": -1}, {"name": "excited", "color": "black", "index": -1}, {"name": "ghz", "color": "black", "index": -1}, {"name": "ground", "color": "black", "index": -1}, {"name": "identity_uncoupled", "color": "black", "index": -1}, {"name": "isdiagonal", "color": "black", "index": -1}, {"name": "jspin", "color": "black", "index": -1}, {"name": "m_degeneracy", "color": "black", "index": -1}, {"name": "num_dicke_ladders", "color": "black", "index": -1}, {"name": "num_dicke_states", "color": "black", "index": -1}, {"name": "num_tls", "color": "black", "index": -1}, {"name": "purity_dicke", "color": "black", "index": -1}, {"name": "spin_algebra", "color": "black", "index": -1}, {"name": "state_degeneracy", "color": "black", "index": -1}, {"name": "superradiant", "color": "black", "index": -1}, {"name": "tau_column", "color": "black", "index": -1}]}, {"name": "solver", "color": "black", "index": -1, "children": [{"name": "ExpectOps", "color": "black", "index": -1}, {"name": "McOptions", "color": "black", "index": -1}, {"name": "Result", "color": "black", "index": -1}, {"name": "SolverConfiguration", "color": "black", "index": -1}, {"name": "SolverOptions", "color": "black", "index": -1}, {"name": "SolverSystem", "color": "black", "index": -1}, {"name": "Stats", "color": "black", "index": -1}]}, {"name": "pdpsolve", "color": "black", "index": -1, "children": [{"name": "StochasticSolverOptions", "color": "black", "index": -1}, {"name": "main_smepdpsolve", "color": "black", "index": -1}, {"name": "main_ssepdpsolve", "color": "black", "index": -1}]}, {"name": "_mcsolve", "color": "black", "index": -1, "children": [{"name": "CyMcOde", "color": "black", "index": -1}, {"name": "CyMcOdeDiag", "color": "black", "index": -1}]}, {"name": "_piqs", "color": "black", "index": -1, "children": [{"name": "Dicke", "color": "black", "index": -1}, {"name": "get_blocks", "color": "black", "index": -1}, {"name": "get_index", "color": "black", "index": -1}, {"name": "j_min", "color": "black", "index": -1}, {"name": "j_vals", "color": "black", "index": -1}, {"name": "jmm1_dictionary", "color": "black", "index": -1}, {"name": "m_vals", "color": "black", "index": -1}]}, {"name": "_steadystate", "color": "black", "index": -1, "children": [{"name": "weighted_bipartite_matching", "color": "black", "index": -1}]}, {"name": "_stochastic", "color": "black", "index": -1, "children": [{"name": "GenericSSolver", "color": "black", "index": -1}, {"name": "PcSMESolver", "color": "black", "index": -1}, {"name": "PcSSESolver", "color": "black", "index": -1}, {"name": "PmSMESolver", "color": "black", "index": -1}, {"name": "SMESolver", "color": "black", "index": -1}, {"name": "SSESolver", "color": "black", "index": -1}, {"name": "StochasticSolver", "color": "black", "index": -1}, {"name": "TaylorNoise", "color": "black", "index": -1}, {"name": "normalize_inplace", "color": "black", "index": -1}]}, {"name": "steadystate", "color": "black", "index": -1, "children": [{"name": "build_preconditioner", "color": "black", "index": -1}, {"name": "pseudo_inverse", "color": "black", "index": -1}, {"name": "steadystate", "color": "black", "index": -1}, {"name": "steadystate_floquet", "color": "black", "index": -1}]}, {"name": "correlation", "color": "black", "index": -1, "children": [{"name": "coherence_function_g1", "color": "black", "index": -1}, {"name": "coherence_function_g2", "color": "black", "index": -1}, {"name": "correlation_2op_1t", "color": "black", "index": -1}, {"name": "correlation_2op_2t", "color": "black", "index": -1}, {"name": "correlation_3op_1t", "color": "black", "index": -1}, {"name": "correlation_3op_2t", "color": "black", "index": -1}, {"name": "spectrum", "color": "black", "index": -1}, {"name": "spectrum_correlation_fft", "color": "black", "index": -1}]}, {"name": "mcsolve", "color": "black", "index": -1, "children": [{"name": "mcsolve", "color": "black", "index": -1}, {"name": "qutip_zvode", "color": "black", "index": -1}]}, {"name": "sesolve", "color": "black", "index": -1, "children": [{"name": "sesolve", "color": "black", "index": -1}]}, {"name": "mesolve", "color": "black", "index": -1, "children": [{"name": "mesolve", "color": "black", "index": -1}]}, {"name": "countstat", "color": "black", "index": -1, "children": [{"name": "countstat_current", "color": "black", "index": -1}, {"name": "countstat_current_noise", "color": "black", "index": -1}]}, {"name": "floquet", "color": "black", "index": -1, "children": [{"name": "floquet_basis_transform", "color": "black", "index": -1}, {"name": "floquet_collapse_operators", "color": "black", "index": -1}, {"name": "floquet_markov_mesolve", "color": "black", "index": -1}, {"name": "floquet_master_equation_rates", "color": "black", "index": -1}, {"name": "floquet_master_equation_steadystate", "color": "black", "index": -1}, {"name": "floquet_master_equation_tensor", "color": "black", "index": -1}, {"name": "floquet_modes", "color": "black", "index": -1}, {"name": "floquet_modes_t", "color": "black", "index": -1}, {"name": "floquet_modes_t_lookup", "color": "black", "index": -1}, {"name": "floquet_modes_table", "color": "black", "index": -1}, {"name": "floquet_state_decomposition", "color": "black", "index": -1}, {"name": "floquet_states", "color": "black", "index": -1}, {"name": "floquet_states_t", "color": "black", "index": -1}, {"name": "floquet_wavefunction", "color": "black", "index": -1}, {"name": "floquet_wavefunction_t", "color": "black", "index": -1}, {"name": "fmmesolve", "color": "black", "index": -1}, {"name": "fsesolve", "color": "black", "index": -1}]}, {"name": "stochastic", "color": "black", "index": -1, "children": [{"name": "StochasticSolverOptions", "color": "black", "index": -1}, {"name": "StochasticSolverOptionsPhoto", "color": "black", "index": -1}, {"name": "general_stochastic", "color": "black", "index": -1}, {"name": "photocurrent_mesolve", "color": "black", "index": -1}, {"name": "photocurrent_sesolve", "color": "black", "index": -1}, {"name": "smepdpsolve", "color": "black", "index": -1}, {"name": "smesolve", "color": "black", "index": -1}, {"name": "ssepdpsolve", "color": "black", "index": -1}, {"name": "ssesolve", "color": "black", "index": -1}, {"name": "stochastic_solvers", "color": "black", "index": -1}]}, {"name": "krylovsolve", "color": "black", "index": -1, "children": [{"name": "krylovsolve", "color": "black", "index": -1}, {"name": "lanczos_algorithm", "color": "black", "index": -1}, {"name": "particular_tlist_or_happy_breakdown", "color": "black", "index": -1}]}, {"name": "nonmarkov", "color": "black", "index": -1}, {"name": "rcsolve", "color": "black", "index": -1, "children": [{"name": "rcsolve", "color": "black", "index": -1}]}]}, {"name": "solver", "color": "#043c6b", "index": 1, "children": [{"name": "brmesolve", "color": "#043c6b", "index": 1, "children": [{"name": "BRSolver", "color": "#043c6b", "index": 1}, {"name": "brmesolve", "color": "#043c6b", "index": 1}]}, {"name": "solver_base", "color": "#043c6b", "index": 1, "children": [{"name": "Solver", "color": "#043c6b", "index": 1}]}, {"name": "integrator", "color": "#043c6b", "index": 1, "children": [{"name": "Integrator", "color": "#043c6b", "index": 1}, {"name": "IntegratorException", "color": "#043c6b", "index": 1}]}, {"name": "ode", "color": "#043c6b", "index": 1, "children": [{"name": "qutip_integrator", "color": "#043c6b", "index": 1, "children": [{"name": "IntegratorDiag", "color": "#043c6b", "index": 1}, {"name": "IntegratorVern", "color": "#043c6b", "index": 1}]}, {"name": "explicit_rk", "color": "#043c6b", "index": 1, "children": [{"name": "Explicit_RungeKutta", "color": "#043c6b", "index": 1}]}, {"name": "scipy_integrator", "color": "#043c6b", "index": 1, "children": [{"name": "IntegratorScipyDop853", "color": "#043c6b", "index": 1}, {"name": "IntegratorScipyZvode", "color": "#043c6b", "index": 1}, {"name": "IntegratorScipylsoda", "color": "#043c6b", "index": 1}]}, {"name": "verner7efficient", "color": "#043c6b", "index": 1}, {"name": "verner9efficient", "color": "#043c6b", "index": 1}]}, {"name": "options", "color": "#043c6b", "index": 1, "children": [{"name": "McOptions", "color": "#043c6b", "index": 1}, {"name": "SolverOdeOptions", "color": "#043c6b", "index": 1}, {"name": "SolverOptions", "color": "#043c6b", "index": 1}, {"name": "SolverResultsOptions", "color": "#043c6b", "index": 1}]}, {"name": "result", "color": "#043c6b", "index": 1, "children": [{"name": "ExpectOp", "color": "#043c6b", "index": 1}, {"name": "MultiTrajResult", "color": "#043c6b", "index": 1}, {"name": "MultiTrajResultAveraged", "color": "#043c6b", "index": 1}, {"name": "Result", "color": "#043c6b", "index": 1}]}, {"name": "propagator", "color": "#043c6b", "index": 1, "children": [{"name": "Propagator", "color": "#043c6b", "index": 1}, {"name": "propagator", "color": "#043c6b", "index": 1}, {"name": "propagator_steadystate", "color": "#043c6b", "index": 1}]}, {"name": "mesolve", "color": "#043c6b", "index": 1, "children": [{"name": "MeSolver", "color": "#043c6b", "index": 1}, {"name": "mesolve", "color": "#043c6b", "index": 1}]}, {"name": "sesolve", "color": "#043c6b", "index": 1, "children": [{"name": "SESolver", "color": "#043c6b", "index": 1}, {"name": "sesolve", "color": "#043c6b", "index": 1}]}, {"name": "scattering", "color": "#043c6b", "index": 1, "children": [{"name": "photon_scattering_amplitude", "color": "#043c6b", "index": 1}, {"name": "scattering_probability", "color": "#043c6b", "index": 1}, {"name": "set_partition", "color": "#043c6b", "index": 1}, {"name": "temporal_basis_vector", "color": "#043c6b", "index": 1}, {"name": "temporal_scattered_state", "color": "#043c6b", "index": 1}]}]}, {"name": "settings", "color": "#043c6b", "index": 1, "children": [{"name": "Settings", "color": "#043c6b", "index": 1}, {"name": "available_cpu_count", "color": "#043c6b", "index": 1}]}, {"name": "bloch", "color": "#3f8fd2", "index": 2, "children": [{"name": "Arrow3D", "color": "#3f8fd2", "index": 2}, {"name": "Bloch", "color": "#3f8fd2", "index": 2}]}, {"name": "bloch3d", "color": "#3f8fd2", "index": 2, "children": [{"name": "Bloch3d", "color": "#3f8fd2", "index": 2}]}, {"name": "distributions", "color": "#3f8fd2", "index": 2, "children": [{"name": "Distribution", "color": "#3f8fd2", "index": 2}, {"name": "HarmonicOscillatorProbabilityFunction", "color": "#3f8fd2", "index": 2}, {"name": "HarmonicOscillatorWaveFunction", "color": "#3f8fd2", "index": 2}, {"name": "QDistribution", "color": "#3f8fd2", "index": 2}, {"name": "TwoModeQuadratureCorrelation", "color": "#3f8fd2", "index": 2}, {"name": "WignerDistribution", "color": "#3f8fd2", "index": 2}]}, {"name": "wigner", "color": "#3f8fd2", "index": 2, "children": [{"name": "QFunc", "color": "#3f8fd2", "index": 2}, {"name": "qfunc", "color": "#3f8fd2", "index": 2}, {"name": "spin_q_function", "color": "#3f8fd2", "index": 2}, {"name": "spin_wigner", "color": "#3f8fd2", "index": 2}, {"name": "wigner", "color": "#3f8fd2", "index": 2}, {"name": "wigner_transform", "color": "#3f8fd2", "index": 2}]}, {"name": "visualization", "color": "#3f8fd2", "index": 2, "children": [{"name": "complex_array_to_rgb", "color": "#3f8fd2", "index": 2}, {"name": "energy_level_diagram", "color": "#3f8fd2", "index": 2}, {"name": "fock_distribution", "color": "#3f8fd2", "index": 2}, {"name": "hinton", "color": "#3f8fd2", "index": 2}, {"name": "matrix_histogram", "color": "#3f8fd2", "index": 2}, {"name": "matrix_histogram_complex", "color": "#3f8fd2", "index": 2}, {"name": "plot_energy_levels", "color": "#3f8fd2", "index": 2}, {"name": "plot_expectation_values", "color": "#3f8fd2", "index": 2}, {"name": "plot_fock_distribution", "color": "#3f8fd2", "index": 2}, {"name": "plot_qubism", "color": "#3f8fd2", "index": 2}, {"name": "plot_schmidt", "color": "#3f8fd2", "index": 2}, {"name": "plot_spin_distribution_2d", "color": "#3f8fd2", "index": 2}, {"name": "plot_spin_distribution_3d", "color": "#3f8fd2", "index": 2}, {"name": "plot_wigner", "color": "#3f8fd2", "index": 2}, {"name": "plot_wigner_fock_distribution", "color": "#3f8fd2", "index": 2}, {"name": "plot_wigner_sphere", "color": "#3f8fd2", "index": 2}, {"name": "sphereplot", "color": "#3f8fd2", "index": 2}, {"name": "wigner_fock_distribution", "color": "#3f8fd2", "index": 2}]}, {"name": "orbital", "color": "#3f8fd2", "index": 2, "children": [{"name": "orbital", "color": "#3f8fd2", "index": 2}]}, {"name": "tomography", "color": "#3f8fd2", "index": 2, "children": [{"name": "qpt", "color": "#3f8fd2", "index": 2}, {"name": "qpt_plot", "color": "#3f8fd2", "index": 2}, {"name": "qpt_plot_combined", "color": "#3f8fd2", "index": 2}]}, {"name": "continuous_variables", "color": "#007143", "index": 4, "children": [{"name": "correlation_matrix", "color": "#007143", "index": 4}, {"name": "correlation_matrix_field", "color": "#007143", "index": 4}, {"name": "correlation_matrix_quadrature", "color": "#007143", "index": 4}, {"name": "covariance_matrix", "color": "#007143", "index": 4}, {"name": "logarithmic_negativity", "color": "#007143", "index": 4}, {"name": "wigner_covariance_matrix", "color": "#007143", "index": 4}]}, {"name": "random_objects", "color": "#007143", "index": 4, "children": [{"name": "rand_dm", "color": "#007143", "index": 4}, {"name": "rand_dm_ginibre", "color": "#007143", "index": 4}, {"name": "rand_dm_hs", "color": "#007143", "index": 4}, {"name": "rand_herm", "color": "#007143", "index": 4}, {"name": "rand_jacobi_rotation", "color": "#007143", "index": 4}, {"name": "rand_ket", "color": "#007143", "index": 4}, {"name": "rand_ket_haar", "color": "#007143", "index": 4}, {"name": "rand_kraus_map", "color": "#007143", "index": 4}, {"name": "rand_stochastic", "color": "#007143", "index": 4}, {"name": "rand_super", "color": "#007143", "index": 4}, {"name": "rand_super_bcsz", "color": "#007143", "index": 4}, {"name": "rand_unitary", "color": "#007143", "index": 4}, {"name": "rand_unitary_haar", "color": "#007143", "index": 4}]}, {"name": "three_level_atom", "color": "#007143", "index": 4, "children": [{"name": "three_level_basis", "color": "#007143", "index": 4}, {"name": "three_level_ops", "color": "#007143", "index": 4}]}, {"name": "entropy", "color": "#ff4500", "index": 6, "children": [{"name": "concurrence", "color": "#ff4500", "index": 6}, {"name": "entangling_power", "color": "#ff4500", "index": 6}, {"name": "entropy_conditional", "color": "#ff4500", "index": 6}, {"name": "entropy_linear", "color": "#ff4500", "index": 6}, {"name": "entropy_mutual", "color": "#ff4500", "index": 6}, {"name": "entropy_relative", "color": "#ff4500", "index": 6}, {"name": "entropy_vn", "color": "#ff4500", "index": 6}, {"name": "negativity", "color": "#ff4500", "index": 6}, {"name": "participation_ratio", "color": "#ff4500", "index": 6}]}, {"name": "partial_transpose", "color": "#692102", "index": 7, "children": [{"name": "partial_transpose", "color": "#692102", "index": 7}]}, {"name": "ui", "color": "#bf5730", "index": 8, "children": [{"name": "progressbar", "color": "#bf5730", "index": 8, "children": [{"name": "BaseProgressBar", "color": "#bf5730", "index": 8}, {"name": "EnhancedTextProgressBar", "color": "#bf5730", "index": 8}, {"name": "TextProgressBar", "color": "#bf5730", "index": 8}, {"name": "TqdmProgressBar", "color": "#bf5730", "index": 8}]}]}, {"name": "parallel", "color": "#bf5730", "index": 8, "children": [{"name": "parallel_map", "color": "#bf5730", "index": 8}, {"name": "parfor", "color": "#bf5730", "index": 8}, {"name": "serial_map", "color": "#bf5730", "index": 8}]}, {"name": "utilities", "color": "#bf5730", "index": 8, "children": [{"name": "clebsch", "color": "#bf5730", "index": 8}, {"name": "convert_GHz_to_mK", "color": "#bf5730", "index": 8}, {"name": "convert_GHz_to_meV", "color": "#bf5730", "index": 8}, {"name": "convert_J_to_meV", "color": "#bf5730", "index": 8}, {"name": "convert_mK_to_GHz", "color": "#bf5730", "index": 8}, {"name": "convert_mK_to_meV", "color": "#bf5730", "index": 8}, {"name": "convert_meV_to_GHz", "color": "#bf5730", "index": 8}, {"name": "convert_meV_to_J", "color": "#bf5730", "index": 8}, {"name": "convert_meV_to_mK", "color": "#bf5730", "index": 8}, {"name": "convert_unit", "color": "#bf5730", "index": 8}, {"name": "n_thermal", "color": "#bf5730", "index": 8}]}, {"name": "about", "color": "#bf5730", "index": 8, "children": [{"name": "about", "color": "#bf5730", "index": 8}]}, {"name": "cite", "color": "#bf5730", "index": 8, "children": [{"name": "cite", "color": "#bf5730", "index": 8}]}, {"name": "fileio", "color": "#bf5730", "index": 8, "children": [{"name": "file_data_read", "color": "#bf5730", "index": 8}, {"name": "file_data_store", "color": "#bf5730", "index": 8}, {"name": "qload", "color": "#bf5730", "index": 8}, {"name": "qsave", "color": "#bf5730", "index": 8}]}, {"name": "simdiag", "color": "#bf5730", "index": 8, "children": [{"name": "simdiag", "color": "#bf5730", "index": 8}]}, {"name": "version", "color": "#bf5730", "index": 8}]}
\ No newline at end of file
+{"name": "QuTiP", "color": "black", "index": -1, "children": [{"name": "core", "color": "black", "index": -1, "children": [{"name": "cy", "color": "black", "index": -1, "children": [{"name": "coefficient", "color": "black", "index": -1, "children": [{"name": "Coefficient", "color": "black", "index": -1}, {"name": "ConjCoefficient", "color": "black", "index": -1}, {"name": "FunctionCoefficient", "color": "black", "index": -1}, {"name": "InterCoefficient", "color": "black", "index": -1}, {"name": "MulCoefficient", "color": "black", "index": -1}, {"name": "NormCoefficient", "color": "black", "index": -1}, {"name": "ShiftCoefficient", "color": "black", "index": -1}, {"name": "StrFunctionCoefficient", "color": "black", "index": -1}, {"name": "SumCoefficient", "color": "black", "index": -1}, {"name": "coefficient_function_parameters", "color": "black", "index": -1}, {"name": "proj", "color": "black", "index": -1}]}, {"name": "_element", "color": "black", "index": -1}, {"name": "qobjevo", "color": "black", "index": -1, "children": [{"name": "QobjEvo", "color": "black", "index": -1}]}]}, {"name": "data", "color": "black", "index": -1, "children": [{"name": "csr", "color": "black", "index": -1, "children": [{"name": "CSR", "color": "black", "index": -1}, {"name": "Sorter", "color": "black", "index": -1}, {"name": "copy_structure", "color": "black", "index": -1}, {"name": "diags", "color": "black", "index": -1}, {"name": "empty", "color": "black", "index": -1}, {"name": "empty_like", "color": "black", "index": -1}, {"name": "fast_from_scipy", "color": "black", "index": -1}, {"name": "from_dense", "color": "black", "index": -1}, {"name": "identity", "color": "black", "index": -1}, {"name": "nnz", "color": "black", "index": -1}, {"name": "sorted", "color": "black", "index": -1}, {"name": "zeros", "color": "black", "index": -1}]}, {"name": "base", "color": "black", "index": -1, "children": [{"name": "Data", "color": "black", "index": -1}, {"name": "EfficiencyWarning", "color": "black", "index": -1}]}, {"name": "dense", "color": "black", "index": -1, "children": [{"name": "Dense", "color": "black", "index": -1}, {"name": "OrderEfficiencyWarning", "color": "black", "index": -1}, {"name": "diags", "color": "black", "index": -1}, {"name": "empty", "color": "black", "index": -1}, {"name": "empty_like", "color": "black", "index": -1}, {"name": "fast_from_numpy", "color": "black", "index": -1}, {"name": "from_csr", "color": "black", "index": -1}, {"name": "identity", "color": "black", "index": -1}, {"name": "zeros", "color": "black", "index": -1}]}, {"name": "dispatch", "color": "black", "index": -1, "children": [{"name": "Dispatcher", "color": "black", "index": -1}]}, {"name": "add", "color": "black", "index": -1, "children": [{"name": "add_csr", "color": "black", "index": -1}, {"name": "add_dense", "color": "black", "index": -1}, {"name": "iadd_dense", "color": "black", "index": -1}, {"name": "sub_csr", "color": "black", "index": -1}, {"name": "sub_dense", "color": "black", "index": -1}]}, {"name": "adjoint", "color": "black", "index": -1, "children": [{"name": "adjoint_csr", "color": "black", "index": -1}, {"name": "adjoint_dense", "color": "black", "index": -1}, {"name": "conj_csr", "color": "black", "index": -1}, {"name": "conj_dense", "color": "black", "index": -1}, {"name": "transpose_csr", "color": "black", "index": -1}, {"name": "transpose_dense", "color": "black", "index": -1}]}, {"name": "reshape", "color": "black", "index": -1, "children": [{"name": "column_stack_csr", "color": "black", "index": -1}, {"name": "column_stack_dense", "color": "black", "index": -1}, {"name": "column_unstack_csr", "color": "black", "index": -1}, {"name": "column_unstack_dense", "color": "black", "index": -1}, {"name": "reshape_csr", "color": "black", "index": -1}, {"name": "reshape_dense", "color": "black", "index": -1}, {"name": "split_columns_csr", "color": "black", "index": -1}, {"name": "split_columns_dense", "color": "black", "index": -1}]}, {"name": "constant", "color": "black", "index": -1}, {"name": "convert", "color": "black", "index": -1}, {"name": "eigen", "color": "black", "index": -1, "children": [{"name": "eigs_csr", "color": "black", "index": -1}, {"name": "eigs_dense", "color": "black", "index": -1}]}, {"name": "expect", "color": "black", "index": -1, "children": [{"name": "expect_csr", "color": "black", "index": -1}, {"name": "expect_csr_dense", "color": "black", "index": -1}, {"name": "expect_dense", "color": "black", "index": -1}, {"name": "expect_super_csr", "color": "black", "index": -1}, {"name": "expect_super_csr_dense", "color": "black", "index": -1}, {"name": "expect_super_dense", "color": "black", "index": -1}]}, {"name": "expm", "color": "black", "index": -1, "children": [{"name": "expm_csr", "color": "black", "index": -1}, {"name": "expm_csr_dense", "color": "black", "index": -1}]}, {"name": "properties", "color": "black", "index": -1, "children": [{"name": "isdiag_csr", "color": "black", "index": -1}, {"name": "isherm_csr", "color": "black", "index": -1}, {"name": "iszero_csr", "color": "black", "index": -1}, {"name": "iszero_dense", "color": "black", "index": -1}]}, {"name": "mul", "color": "black", "index": -1, "children": [{"name": "imul_csr", "color": "black", "index": -1}, {"name": "imul_data", "color": "black", "index": -1}, {"name": "imul_dense", "color": "black", "index": -1}, {"name": "mul_csr", "color": "black", "index": -1}, {"name": "mul_dense", "color": "black", "index": -1}, {"name": "neg_csr", "color": "black", "index": -1}, {"name": "neg_dense", "color": "black", "index": -1}]}, {"name": "inner", "color": "black", "index": -1, "children": [{"name": "inner_csr", "color": "black", "index": -1}, {"name": "inner_op_csr", "color": "black", "index": -1}]}, {"name": "linalg", "color": "black", "index": -1, "children": [{"name": "inv_csr", "color": "black", "index": -1}, {"name": "inv_dense", "color": "black", "index": -1}]}, {"name": "kron", "color": "black", "index": -1, "children": [{"name": "kron_csr", "color": "black", "index": -1}, {"name": "kron_dense", "color": "black", "index": -1}]}, {"name": "make", "color": "black", "index": -1, "children": [{"name": "one_element_csr", "color": "black", "index": -1}, {"name": "one_element_dense", "color": "black", "index": -1}]}, {"name": "matmul", "color": "black", "index": -1, "children": [{"name": "matmul_csr", "color": "black", "index": -1}, {"name": "matmul_csr_dense_dense", "color": "black", "index": -1}, {"name": "matmul_dense", "color": "black", "index": -1}, {"name": "multiply_csr", "color": "black", "index": -1}, {"name": "multiply_dense", "color": "black", "index": -1}]}, {"name": "norm", "color": "black", "index": -1, "children": [{"name": "frobenius_csr", "color": "black", "index": -1}, {"name": "frobenius_data", "color": "black", "index": -1}, {"name": "frobenius_dense", "color": "black", "index": -1}, {"name": "l2_csr", "color": "black", "index": -1}, {"name": "l2_dense", "color": "black", "index": -1}, {"name": "max_csr", "color": "black", "index": -1}, {"name": "max_dense", "color": "black", "index": -1}, {"name": "one_csr", "color": "black", "index": -1}, {"name": "one_dense", "color": "black", "index": -1}, {"name": "trace_csr", "color": "black", "index": -1}, {"name": "trace_dense", "color": "black", "index": -1}]}, {"name": "permute", "color": "black", "index": -1, "children": [{"name": "dimensions_csr", "color": "black", "index": -1}, {"name": "indices_csr", "color": "black", "index": -1}]}, {"name": "pow", "color": "black", "index": -1, "children": [{"name": "pow_csr", "color": "black", "index": -1}]}, {"name": "project", "color": "black", "index": -1, "children": [{"name": "project_csr", "color": "black", "index": -1}, {"name": "project_dense", "color": "black", "index": -1}]}, {"name": "ptrace", "color": "black", "index": -1, "children": [{"name": "ptrace_csr", "color": "black", "index": -1}, {"name": "ptrace_csr_dense", "color": "black", "index": -1}, {"name": "ptrace_dense", "color": "black", "index": -1}]}, {"name": "tidyup", "color": "black", "index": -1, "children": [{"name": "tidyup_csr", "color": "black", "index": -1}, {"name": "tidyup_dense", "color": "black", "index": -1}]}, {"name": "trace", "color": "black", "index": -1, "children": [{"name": "trace_csr", "color": "black", "index": -1}, {"name": "trace_dense", "color": "black", "index": -1}]}]}, {"name": "coefficient", "color": "black", "index": -1, "children": [{"name": "CompilationOptions", "color": "black", "index": -1}, {"name": "StringParsingWarning", "color": "black", "index": -1}, {"name": "clean_compiled_coefficient", "color": "black", "index": -1}, {"name": "coeff_from_str", "color": "black", "index": -1}, {"name": "coefficient", "color": "black", "index": -1}, {"name": "compileType", "color": "black", "index": -1}, {"name": "compile_code", "color": "black", "index": -1}, {"name": "conj", "color": "black", "index": -1}, {"name": "extract_constant", "color": "black", "index": -1}, {"name": "extract_cte_pattern", "color": "black", "index": -1}, {"name": "find_type_from_str", "color": "black", "index": -1}, {"name": "fix_type", "color": "black", "index": -1}, {"name": "fromstr", "color": "black", "index": -1}, {"name": "make_cy_code", "color": "black", "index": -1}, {"name": "norm", "color": "black", "index": -1}, {"name": "parse", "color": "black", "index": -1}, {"name": "proj", "color": "black", "index": -1}, {"name": "shift", "color": "black", "index": -1}, {"name": "space_parts", "color": "black", "index": -1}, {"name": "test_parsed", "color": "black", "index": -1}, {"name": "try_import", "color": "black", "index": -1}, {"name": "try_parse", "color": "black", "index": -1}, {"name": "use_hinted_type", "color": "black", "index": -1}]}, {"name": "options", "color": "black", "index": -1, "children": [{"name": "CoreOptions", "color": "black", "index": -1}, {"name": "QutipOptions", "color": "black", "index": -1}]}, {"name": "qobj", "color": "black", "index": -1, "children": [{"name": "Qobj", "color": "black", "index": -1}, {"name": "isbra", "color": "black", "index": -1}, {"name": "isherm", "color": "black", "index": -1}, {"name": "isket", "color": "black", "index": -1}, {"name": "isoper", "color": "black", "index": -1}, {"name": "isoperbra", "color": "black", "index": -1}, {"name": "isoperket", "color": "black", "index": -1}, {"name": "issuper", "color": "black", "index": -1}, {"name": "ptrace", "color": "black", "index": -1}]}, {"name": "dimensions", "color": "black", "index": -1, "children": [{"name": "collapse_dims_oper", "color": "black", "index": -1}, {"name": "collapse_dims_super", "color": "black", "index": -1}, {"name": "deep_map", "color": "black", "index": -1}, {"name": "deep_remove", "color": "black", "index": -1}, {"name": "dims_idxs_to_tensor_idxs", "color": "black", "index": -1}, {"name": "dims_to_tensor_perm", "color": "black", "index": -1}, {"name": "dims_to_tensor_shape", "color": "black", "index": -1}, {"name": "enumerate_flat", "color": "black", "index": -1}, {"name": "flatten", "color": "black", "index": -1}, {"name": "is_scalar", "color": "black", "index": -1}, {"name": "is_vector", "color": "black", "index": -1}, {"name": "is_vectorized_oper", "color": "black", "index": -1}, {"name": "type_from_dims", "color": "black", "index": -1}, {"name": "unflatten", "color": "black", "index": -1}]}, {"name": "metrics", "color": "black", "index": -1, "children": [{"name": "average_gate_fidelity", "color": "black", "index": -1}, {"name": "bures_angle", "color": "black", "index": -1}, {"name": "bures_dist", "color": "black", "index": -1}, {"name": "dnorm", "color": "black", "index": -1}, {"name": "fidelity", "color": "black", "index": -1}, {"name": "hellinger_dist", "color": "black", "index": -1}, {"name": "hilbert_dist", "color": "black", "index": -1}, {"name": "process_fidelity", "color": "black", "index": -1}, {"name": "tracedist", "color": "black", "index": -1}, {"name": "unitarity", "color": "black", "index": -1}]}, {"name": "superop_reps", "color": "black", "index": -1, "children": [{"name": "isqubitdims", "color": "black", "index": -1}, {"name": "kraus_to_choi", "color": "black", "index": -1}, {"name": "kraus_to_super", "color": "black", "index": -1}, {"name": "to_chi", "color": "black", "index": -1}, {"name": "to_choi", "color": "black", "index": -1}, {"name": "to_kraus", "color": "black", "index": -1}, {"name": "to_stinespring", "color": "black", "index": -1}, {"name": "to_super", "color": "black", "index": -1}]}, {"name": "states", "color": "black", "index": -1, "children": [{"name": "basis", "color": "black", "index": -1}, {"name": "bell_state", "color": "black", "index": -1}, {"name": "bra", "color": "black", "index": -1}, {"name": "coherent", "color": "black", "index": -1}, {"name": "coherent_dm", "color": "black", "index": -1}, {"name": "enr_fock", "color": "black", "index": -1}, {"name": "enr_state_dictionaries", "color": "black", "index": -1}, {"name": "enr_thermal_dm", "color": "black", "index": -1}, {"name": "fock", "color": "black", "index": -1}, {"name": "fock_dm", "color": "black", "index": -1}, {"name": "ghz_state", "color": "black", "index": -1}, {"name": "ket", "color": "black", "index": -1}, {"name": "ket2dm", "color": "black", "index": -1}, {"name": "maximally_mixed_dm", "color": "black", "index": -1}, {"name": "phase_basis", "color": "black", "index": -1}, {"name": "projection", "color": "black", "index": -1}, {"name": "qstate", "color": "black", "index": -1}, {"name": "qutrit_basis", "color": "black", "index": -1}, {"name": "singlet_state", "color": "black", "index": -1}, {"name": "spin_coherent", "color": "black", "index": -1}, {"name": "spin_state", "color": "black", "index": -1}, {"name": "state_index_number", "color": "black", "index": -1}, {"name": "state_number_enumerate", "color": "black", "index": -1}, {"name": "state_number_index", "color": "black", "index": -1}, {"name": "state_number_qobj", "color": "black", "index": -1}, {"name": "thermal_dm", "color": "black", "index": -1}, {"name": "triplet_states", "color": "black", "index": -1}, {"name": "w_state", "color": "black", "index": -1}, {"name": "zero_ket", "color": "black", "index": -1}]}, {"name": "operators", "color": "black", "index": -1, "children": [{"name": "charge", "color": "black", "index": -1}, {"name": "commutator", "color": "black", "index": -1}, {"name": "create", "color": "black", "index": -1}, {"name": "destroy", "color": "black", "index": -1}, {"name": "displace", "color": "black", "index": -1}, {"name": "enr_destroy", "color": "black", "index": -1}, {"name": "enr_identity", "color": "black", "index": -1}, {"name": "qeye", "color": "black", "index": -1}, {"name": "jmat", "color": "black", "index": -1}, {"name": "momentum", "color": "black", "index": -1}, {"name": "num", "color": "black", "index": -1}, {"name": "phase", "color": "black", "index": -1}, {"name": "position", "color": "black", "index": -1}, {"name": "qdiags", "color": "black", "index": -1}, {"name": "qft", "color": "black", "index": -1}, {"name": "qutrit_ops", "color": "black", "index": -1}, {"name": "qzero", "color": "black", "index": -1}, {"name": "sigmam", "color": "black", "index": -1}, {"name": "sigmap", "color": "black", "index": -1}, {"name": "sigmax", "color": "black", "index": -1}, {"name": "sigmay", "color": "black", "index": -1}, {"name": "sigmaz", "color": "black", "index": -1}, {"name": "spin_J_set", "color": "black", "index": -1}, {"name": "spin_Jm", "color": "black", "index": -1}, {"name": "spin_Jp", "color": "black", "index": -1}, {"name": "spin_Jx", "color": "black", "index": -1}, {"name": "spin_Jy", "color": "black", "index": -1}, {"name": "spin_Jz", "color": "black", "index": -1}, {"name": "squeeze", "color": "black", "index": -1}, {"name": "squeezing", "color": "black", "index": -1}, {"name": "tunneling", "color": "black", "index": -1}]}, {"name": "tensor", "color": "black", "index": -1, "children": [{"name": "composite", "color": "black", "index": -1}, {"name": "expand_operator", "color": "black", "index": -1}, {"name": "super_tensor", "color": "black", "index": -1}, {"name": "tensor", "color": "black", "index": -1}, {"name": "tensor_contract", "color": "black", "index": -1}, {"name": "tensor_swap", "color": "black", "index": -1}]}, {"name": "superoperator", "color": "black", "index": -1, "children": [{"name": "lindblad_dissipator", "color": "black", "index": -1}, {"name": "liouvillian", "color": "black", "index": -1}, {"name": "operator_to_vector", "color": "black", "index": -1}, {"name": "reshuffle", "color": "black", "index": -1}, {"name": "spost", "color": "black", "index": -1}, {"name": "spre", "color": "black", "index": -1}, {"name": "sprepost", "color": "black", "index": -1}, {"name": "stack_columns", "color": "black", "index": -1}, {"name": "stacked_index", "color": "black", "index": -1}, {"name": "unstack_columns", "color": "black", "index": -1}, {"name": "unstacked_index", "color": "black", "index": -1}, {"name": "vector_to_operator", "color": "black", "index": -1}]}, {"name": "semidefinite", "color": "black", "index": -1, "children": [{"name": "Complex", "color": "black", "index": -1}, {"name": "bmat", "color": "black", "index": -1}, {"name": "complex_var", "color": "black", "index": -1}, {"name": "conj", "color": "black", "index": -1}, {"name": "dag", "color": "black", "index": -1}, {"name": "dens", "color": "black", "index": -1}, {"name": "dnorm_problem", "color": "black", "index": -1}, {"name": "dnorm_sparse_problem", "color": "black", "index": -1}, {"name": "herm", "color": "black", "index": -1}, {"name": "initialize_constraints_on_dnorm_problem", "color": "black", "index": -1}, {"name": "kron", "color": "black", "index": -1}, {"name": "memoize", "color": "black", "index": -1}, {"name": "pos", "color": "black", "index": -1}, {"name": "pos_noherm", "color": "black", "index": -1}, {"name": "qudit_swap", "color": "black", "index": -1}]}, {"name": "_brtensor", "color": "black", "index": -1}, {"name": "_brtools", "color": "black", "index": -1, "children": [{"name": "SpectraCoefficient", "color": "black", "index": -1}, {"name": "matmul_var_data", "color": "black", "index": -1}]}, {"name": "blochredfield", "color": "black", "index": -1, "children": [{"name": "bloch_redfield_tensor", "color": "black", "index": -1}, {"name": "brterm", "color": "black", "index": -1}]}, {"name": "expect", "color": "black", "index": -1, "children": [{"name": "expect", "color": "black", "index": -1}, {"name": "variance", "color": "black", "index": -1}]}, {"name": "gates", "color": "black", "index": -1, "children": [{"name": "berkeley", "color": "black", "index": -1}, {"name": "cnot", "color": "black", "index": -1}, {"name": "cphase", "color": "black", "index": -1}, {"name": "cs_gate", "color": "black", "index": -1}, {"name": "csign", "color": "black", "index": -1}, {"name": "ct_gate", "color": "black", "index": -1}, {"name": "cy_gate", "color": "black", "index": -1}, {"name": "cz_gate", "color": "black", "index": -1}, {"name": "fredkin", "color": "black", "index": -1}, {"name": "globalphase", "color": "black", "index": -1}, {"name": "hadamard_transform", "color": "black", "index": -1}, {"name": "iswap", "color": "black", "index": -1}, {"name": "molmer_sorensen", "color": "black", "index": -1}, {"name": "phasegate", "color": "black", "index": -1}, {"name": "qrot", "color": "black", "index": -1}, {"name": "qubit_clifford_group", "color": "black", "index": -1}, {"name": "rx", "color": "black", "index": -1}, {"name": "ry", "color": "black", "index": -1}, {"name": "rz", "color": "black", "index": -1}, {"name": "s_gate", "color": "black", "index": -1}, {"name": "snot", "color": "black", "index": -1}, {"name": "sqrtiswap", "color": "black", "index": -1}, {"name": "sqrtnot", "color": "black", "index": -1}, {"name": "sqrtswap", "color": "black", "index": -1}, {"name": "swap", "color": "black", "index": -1}, {"name": "swapalpha", "color": "black", "index": -1}, {"name": "t_gate", "color": "black", "index": -1}, {"name": "toffoli", "color": "black", "index": -1}]}, {"name": "subsystem_apply", "color": "black", "index": -1, "children": [{"name": "subsystem_apply", "color": "black", "index": -1}]}]}, {"name": "optionsclass", "color": "black", "index": -1, "children": [{"name": "AllOptions", "color": "black", "index": -1}, {"name": "optionsclass", "color": "black", "index": -1}]}, {"name": "solve", "color": "black", "index": -1, "children": [{"name": "piqs", "color": "black", "index": -1, "children": [{"name": "Dicke", "color": "black", "index": -1}, {"name": "Pim", "color": "black", "index": -1}, {"name": "am", "color": "black", "index": -1}, {"name": "ap", "color": "black", "index": -1}, {"name": "block_matrix", "color": "black", "index": -1}, {"name": "collapse_uncoupled", "color": "black", "index": -1}, {"name": "css", "color": "black", "index": -1}, {"name": "dicke", "color": "black", "index": -1}, {"name": "dicke_basis", "color": "black", "index": -1}, {"name": "dicke_blocks", "color": "black", "index": -1}, {"name": "dicke_blocks_full", "color": "black", "index": -1}, {"name": "dicke_function_trace", "color": "black", "index": -1}, {"name": "energy_degeneracy", "color": "black", "index": -1}, {"name": "entropy_vn_dicke", "color": "black", "index": -1}, {"name": "excited", "color": "black", "index": -1}, {"name": "ghz", "color": "black", "index": -1}, {"name": "ground", "color": "black", "index": -1}, {"name": "identity_uncoupled", "color": "black", "index": -1}, {"name": "isdiagonal", "color": "black", "index": -1}, {"name": "jspin", "color": "black", "index": -1}, {"name": "m_degeneracy", "color": "black", "index": -1}, {"name": "num_dicke_ladders", "color": "black", "index": -1}, {"name": "num_dicke_states", "color": "black", "index": -1}, {"name": "num_tls", "color": "black", "index": -1}, {"name": "purity_dicke", "color": "black", "index": -1}, {"name": "spin_algebra", "color": "black", "index": -1}, {"name": "state_degeneracy", "color": "black", "index": -1}, {"name": "superradiant", "color": "black", "index": -1}, {"name": "tau_column", "color": "black", "index": -1}]}, {"name": "solver", "color": "black", "index": -1, "children": [{"name": "ExpectOps", "color": "black", "index": -1}, {"name": "McOptions", "color": "black", "index": -1}, {"name": "Result", "color": "black", "index": -1}, {"name": "SolverConfiguration", "color": "black", "index": -1}, {"name": "SolverOptions", "color": "black", "index": -1}, {"name": "SolverSystem", "color": "black", "index": -1}, {"name": "Stats", "color": "black", "index": -1}]}, {"name": "pdpsolve", "color": "black", "index": -1, "children": [{"name": "StochasticSolverOptions", "color": "black", "index": -1}, {"name": "main_smepdpsolve", "color": "black", "index": -1}, {"name": "main_ssepdpsolve", "color": "black", "index": -1}]}, {"name": "_mcsolve", "color": "black", "index": -1, "children": [{"name": "CyMcOde", "color": "black", "index": -1}, {"name": "CyMcOdeDiag", "color": "black", "index": -1}]}, {"name": "_piqs", "color": "black", "index": -1, "children": [{"name": "Dicke", "color": "black", "index": -1}, {"name": "get_blocks", "color": "black", "index": -1}, {"name": "get_index", "color": "black", "index": -1}, {"name": "j_min", "color": "black", "index": -1}, {"name": "j_vals", "color": "black", "index": -1}, {"name": "jmm1_dictionary", "color": "black", "index": -1}, {"name": "m_vals", "color": "black", "index": -1}]}, {"name": "_steadystate", "color": "black", "index": -1, "children": [{"name": "weighted_bipartite_matching", "color": "black", "index": -1}]}, {"name": "_stochastic", "color": "black", "index": -1, "children": [{"name": "GenericSSolver", "color": "black", "index": -1}, {"name": "PcSMESolver", "color": "black", "index": -1}, {"name": "PcSSESolver", "color": "black", "index": -1}, {"name": "PmSMESolver", "color": "black", "index": -1}, {"name": "SMESolver", "color": "black", "index": -1}, {"name": "SSESolver", "color": "black", "index": -1}, {"name": "StochasticSolver", "color": "black", "index": -1}, {"name": "TaylorNoise", "color": "black", "index": -1}, {"name": "normalize_inplace", "color": "black", "index": -1}]}, {"name": "steadystate", "color": "black", "index": -1, "children": [{"name": "build_preconditioner", "color": "black", "index": -1}, {"name": "pseudo_inverse", "color": "black", "index": -1}, {"name": "steadystate", "color": "black", "index": -1}, {"name": "steadystate_floquet", "color": "black", "index": -1}]}, {"name": "correlation", "color": "black", "index": -1, "children": [{"name": "coherence_function_g1", "color": "black", "index": -1}, {"name": "coherence_function_g2", "color": "black", "index": -1}, {"name": "correlation_2op_1t", "color": "black", "index": -1}, {"name": "correlation_2op_2t", "color": "black", "index": -1}, {"name": "correlation_3op_1t", "color": "black", "index": -1}, {"name": "correlation_3op_2t", "color": "black", "index": -1}, {"name": "spectrum", "color": "black", "index": -1}, {"name": "spectrum_correlation_fft", "color": "black", "index": -1}]}, {"name": "mcsolve", "color": "black", "index": -1, "children": [{"name": "mcsolve", "color": "black", "index": -1}, {"name": "qutip_zvode", "color": "black", "index": -1}]}, {"name": "sesolve", "color": "black", "index": -1, "children": [{"name": "sesolve", "color": "black", "index": -1}]}, {"name": "mesolve", "color": "black", "index": -1, "children": [{"name": "mesolve", "color": "black", "index": -1}]}, {"name": "countstat", "color": "black", "index": -1, "children": [{"name": "countstat_current", "color": "black", "index": -1}, {"name": "countstat_current_noise", "color": "black", "index": -1}]}, {"name": "floquet", "color": "black", "index": -1, "children": [{"name": "floquet_basis_transform", "color": "black", "index": -1}, {"name": "floquet_collapse_operators", "color": "black", "index": -1}, {"name": "floquet_markov_mesolve", "color": "black", "index": -1}, {"name": "floquet_master_equation_rates", "color": "black", "index": -1}, {"name": "floquet_master_equation_steadystate", "color": "black", "index": -1}, {"name": "floquet_master_equation_tensor", "color": "black", "index": -1}, {"name": "floquet_modes", "color": "black", "index": -1}, {"name": "floquet_modes_t", "color": "black", "index": -1}, {"name": "floquet_modes_t_lookup", "color": "black", "index": -1}, {"name": "floquet_modes_table", "color": "black", "index": -1}, {"name": "floquet_state_decomposition", "color": "black", "index": -1}, {"name": "floquet_states", "color": "black", "index": -1}, {"name": "floquet_states_t", "color": "black", "index": -1}, {"name": "floquet_wavefunction", "color": "black", "index": -1}, {"name": "floquet_wavefunction_t", "color": "black", "index": -1}, {"name": "fmmesolve", "color": "black", "index": -1}, {"name": "fsesolve", "color": "black", "index": -1}]}, {"name": "stochastic", "color": "black", "index": -1, "children": [{"name": "StochasticSolverOptions", "color": "black", "index": -1}, {"name": "StochasticSolverOptionsPhoto", "color": "black", "index": -1}, {"name": "general_stochastic", "color": "black", "index": -1}, {"name": "photocurrent_mesolve", "color": "black", "index": -1}, {"name": "photocurrent_sesolve", "color": "black", "index": -1}, {"name": "smepdpsolve", "color": "black", "index": -1}, {"name": "smesolve", "color": "black", "index": -1}, {"name": "ssepdpsolve", "color": "black", "index": -1}, {"name": "ssesolve", "color": "black", "index": -1}, {"name": "stochastic_solvers", "color": "black", "index": -1}]}, {"name": "krylovsolve", "color": "black", "index": -1, "children": [{"name": "krylovsolve", "color": "black", "index": -1}, {"name": "lanczos_algorithm", "color": "black", "index": -1}, {"name": "particular_tlist_or_happy_breakdown", "color": "black", "index": -1}]}, {"name": "nonmarkov", "color": "black", "index": -1}, {"name": "rcsolve", "color": "black", "index": -1, "children": [{"name": "rcsolve", "color": "black", "index": -1}]}]}, {"name": "solver", "color": "#043c6b", "index": 1, "children": [{"name": "brmesolve", "color": "#043c6b", "index": 1, "children": [{"name": "BRSolver", "color": "#043c6b", "index": 1}, {"name": "brmesolve", "color": "#043c6b", "index": 1}]}, {"name": "solver_base", "color": "#043c6b", "index": 1, "children": [{"name": "Solver", "color": "#043c6b", "index": 1}]}, {"name": "integrator", "color": "#043c6b", "index": 1, "children": [{"name": "Integrator", "color": "#043c6b", "index": 1}, {"name": "IntegratorException", "color": "#043c6b", "index": 1}]}, {"name": "ode", "color": "#043c6b", "index": 1, "children": [{"name": "qutip_integrator", "color": "#043c6b", "index": 1, "children": [{"name": "IntegratorDiag", "color": "#043c6b", "index": 1}, {"name": "IntegratorVern", "color": "#043c6b", "index": 1}]}, {"name": "explicit_rk", "color": "#043c6b", "index": 1, "children": [{"name": "Explicit_RungeKutta", "color": "#043c6b", "index": 1}]}, {"name": "scipy_integrator", "color": "#043c6b", "index": 1, "children": [{"name": "IntegratorScipyDop853", "color": "#043c6b", "index": 1}, {"name": "IntegratorScipyZvode", "color": "#043c6b", "index": 1}, {"name": "IntegratorScipylsoda", "color": "#043c6b", "index": 1}]}, {"name": "verner7efficient", "color": "#043c6b", "index": 1}, {"name": "verner9efficient", "color": "#043c6b", "index": 1}]}, {"name": "options", "color": "#043c6b", "index": 1, "children": [{"name": "McOptions", "color": "#043c6b", "index": 1}, {"name": "SolverOdeOptions", "color": "#043c6b", "index": 1}, {"name": "SolverOptions", "color": "#043c6b", "index": 1}, {"name": "SolverResultsOptions", "color": "#043c6b", "index": 1}]}, {"name": "result", "color": "#043c6b", "index": 1, "children": [{"name": "ExpectOp", "color": "#043c6b", "index": 1}, {"name": "MultiTrajResult", "color": "#043c6b", "index": 1}, {"name": "MultiTrajResultAveraged", "color": "#043c6b", "index": 1}, {"name": "Result", "color": "#043c6b", "index": 1}]}, {"name": "propagator", "color": "#043c6b", "index": 1, "children": [{"name": "Propagator", "color": "#043c6b", "index": 1}, {"name": "propagator", "color": "#043c6b", "index": 1}, {"name": "propagator_steadystate", "color": "#043c6b", "index": 1}]}, {"name": "mesolve", "color": "#043c6b", "index": 1, "children": [{"name": "MESolver", "color": "#043c6b", "index": 1}, {"name": "mesolve", "color": "#043c6b", "index": 1}]}, {"name": "sesolve", "color": "#043c6b", "index": 1, "children": [{"name": "SESolver", "color": "#043c6b", "index": 1}, {"name": "sesolve", "color": "#043c6b", "index": 1}]}, {"name": "scattering", "color": "#043c6b", "index": 1, "children": [{"name": "photon_scattering_amplitude", "color": "#043c6b", "index": 1}, {"name": "scattering_probability", "color": "#043c6b", "index": 1}, {"name": "set_partition", "color": "#043c6b", "index": 1}, {"name": "temporal_basis_vector", "color": "#043c6b", "index": 1}, {"name": "temporal_scattered_state", "color": "#043c6b", "index": 1}]}]}, {"name": "settings", "color": "#043c6b", "index": 1, "children": [{"name": "Settings", "color": "#043c6b", "index": 1}, {"name": "available_cpu_count", "color": "#043c6b", "index": 1}]}, {"name": "bloch", "color": "#3f8fd2", "index": 2, "children": [{"name": "Arrow3D", "color": "#3f8fd2", "index": 2}, {"name": "Bloch", "color": "#3f8fd2", "index": 2}]}, {"name": "bloch3d", "color": "#3f8fd2", "index": 2, "children": [{"name": "Bloch3d", "color": "#3f8fd2", "index": 2}]}, {"name": "distributions", "color": "#3f8fd2", "index": 2, "children": [{"name": "Distribution", "color": "#3f8fd2", "index": 2}, {"name": "HarmonicOscillatorProbabilityFunction", "color": "#3f8fd2", "index": 2}, {"name": "HarmonicOscillatorWaveFunction", "color": "#3f8fd2", "index": 2}, {"name": "QDistribution", "color": "#3f8fd2", "index": 2}, {"name": "TwoModeQuadratureCorrelation", "color": "#3f8fd2", "index": 2}, {"name": "WignerDistribution", "color": "#3f8fd2", "index": 2}]}, {"name": "wigner", "color": "#3f8fd2", "index": 2, "children": [{"name": "QFunc", "color": "#3f8fd2", "index": 2}, {"name": "qfunc", "color": "#3f8fd2", "index": 2}, {"name": "spin_q_function", "color": "#3f8fd2", "index": 2}, {"name": "spin_wigner", "color": "#3f8fd2", "index": 2}, {"name": "wigner", "color": "#3f8fd2", "index": 2}, {"name": "wigner_transform", "color": "#3f8fd2", "index": 2}]}, {"name": "visualization", "color": "#3f8fd2", "index": 2, "children": [{"name": "complex_array_to_rgb", "color": "#3f8fd2", "index": 2}, {"name": "energy_level_diagram", "color": "#3f8fd2", "index": 2}, {"name": "fock_distribution", "color": "#3f8fd2", "index": 2}, {"name": "hinton", "color": "#3f8fd2", "index": 2}, {"name": "matrix_histogram", "color": "#3f8fd2", "index": 2}, {"name": "matrix_histogram_complex", "color": "#3f8fd2", "index": 2}, {"name": "plot_energy_levels", "color": "#3f8fd2", "index": 2}, {"name": "plot_expectation_values", "color": "#3f8fd2", "index": 2}, {"name": "plot_fock_distribution", "color": "#3f8fd2", "index": 2}, {"name": "plot_qubism", "color": "#3f8fd2", "index": 2}, {"name": "plot_schmidt", "color": "#3f8fd2", "index": 2}, {"name": "plot_spin_distribution_2d", "color": "#3f8fd2", "index": 2}, {"name": "plot_spin_distribution_3d", "color": "#3f8fd2", "index": 2}, {"name": "plot_wigner", "color": "#3f8fd2", "index": 2}, {"name": "plot_wigner_fock_distribution", "color": "#3f8fd2", "index": 2}, {"name": "plot_wigner_sphere", "color": "#3f8fd2", "index": 2}, {"name": "sphereplot", "color": "#3f8fd2", "index": 2}, {"name": "wigner_fock_distribution", "color": "#3f8fd2", "index": 2}]}, {"name": "orbital", "color": "#3f8fd2", "index": 2, "children": [{"name": "orbital", "color": "#3f8fd2", "index": 2}]}, {"name": "tomography", "color": "#3f8fd2", "index": 2, "children": [{"name": "qpt", "color": "#3f8fd2", "index": 2}, {"name": "qpt_plot", "color": "#3f8fd2", "index": 2}, {"name": "qpt_plot_combined", "color": "#3f8fd2", "index": 2}]}, {"name": "continuous_variables", "color": "#007143", "index": 4, "children": [{"name": "correlation_matrix", "color": "#007143", "index": 4}, {"name": "correlation_matrix_field", "color": "#007143", "index": 4}, {"name": "correlation_matrix_quadrature", "color": "#007143", "index": 4}, {"name": "covariance_matrix", "color": "#007143", "index": 4}, {"name": "logarithmic_negativity", "color": "#007143", "index": 4}, {"name": "wigner_covariance_matrix", "color": "#007143", "index": 4}]}, {"name": "random_objects", "color": "#007143", "index": 4, "children": [{"name": "rand_dm", "color": "#007143", "index": 4}, {"name": "rand_dm_ginibre", "color": "#007143", "index": 4}, {"name": "rand_dm_hs", "color": "#007143", "index": 4}, {"name": "rand_herm", "color": "#007143", "index": 4}, {"name": "rand_jacobi_rotation", "color": "#007143", "index": 4}, {"name": "rand_ket", "color": "#007143", "index": 4}, {"name": "rand_ket_haar", "color": "#007143", "index": 4}, {"name": "rand_kraus_map", "color": "#007143", "index": 4}, {"name": "rand_stochastic", "color": "#007143", "index": 4}, {"name": "rand_super", "color": "#007143", "index": 4}, {"name": "rand_super_bcsz", "color": "#007143", "index": 4}, {"name": "rand_unitary", "color": "#007143", "index": 4}, {"name": "rand_unitary_haar", "color": "#007143", "index": 4}]}, {"name": "three_level_atom", "color": "#007143", "index": 4, "children": [{"name": "three_level_basis", "color": "#007143", "index": 4}, {"name": "three_level_ops", "color": "#007143", "index": 4}]}, {"name": "entropy", "color": "#ff4500", "index": 6, "children": [{"name": "concurrence", "color": "#ff4500", "index": 6}, {"name": "entangling_power", "color": "#ff4500", "index": 6}, {"name": "entropy_conditional", "color": "#ff4500", "index": 6}, {"name": "entropy_linear", "color": "#ff4500", "index": 6}, {"name": "entropy_mutual", "color": "#ff4500", "index": 6}, {"name": "entropy_relative", "color": "#ff4500", "index": 6}, {"name": "entropy_vn", "color": "#ff4500", "index": 6}, {"name": "negativity", "color": "#ff4500", "index": 6}, {"name": "participation_ratio", "color": "#ff4500", "index": 6}]}, {"name": "partial_transpose", "color": "#692102", "index": 7, "children": [{"name": "partial_transpose", "color": "#692102", "index": 7}]}, {"name": "ui", "color": "#bf5730", "index": 8, "children": [{"name": "progressbar", "color": "#bf5730", "index": 8, "children": [{"name": "BaseProgressBar", "color": "#bf5730", "index": 8}, {"name": "EnhancedTextProgressBar", "color": "#bf5730", "index": 8}, {"name": "TextProgressBar", "color": "#bf5730", "index": 8}, {"name": "TqdmProgressBar", "color": "#bf5730", "index": 8}]}]}, {"name": "parallel", "color": "#bf5730", "index": 8, "children": [{"name": "parallel_map", "color": "#bf5730", "index": 8}, {"name": "parfor", "color": "#bf5730", "index": 8}, {"name": "serial_map", "color": "#bf5730", "index": 8}]}, {"name": "utilities", "color": "#bf5730", "index": 8, "children": [{"name": "clebsch", "color": "#bf5730", "index": 8}, {"name": "convert_GHz_to_mK", "color": "#bf5730", "index": 8}, {"name": "convert_GHz_to_meV", "color": "#bf5730", "index": 8}, {"name": "convert_J_to_meV", "color": "#bf5730", "index": 8}, {"name": "convert_mK_to_GHz", "color": "#bf5730", "index": 8}, {"name": "convert_mK_to_meV", "color": "#bf5730", "index": 8}, {"name": "convert_meV_to_GHz", "color": "#bf5730", "index": 8}, {"name": "convert_meV_to_J", "color": "#bf5730", "index": 8}, {"name": "convert_meV_to_mK", "color": "#bf5730", "index": 8}, {"name": "convert_unit", "color": "#bf5730", "index": 8}, {"name": "n_thermal", "color": "#bf5730", "index": 8}]}, {"name": "about", "color": "#bf5730", "index": 8, "children": [{"name": "about", "color": "#bf5730", "index": 8}]}, {"name": "cite", "color": "#bf5730", "index": 8, "children": [{"name": "cite", "color": "#bf5730", "index": 8}]}, {"name": "fileio", "color": "#bf5730", "index": 8, "children": [{"name": "file_data_read", "color": "#bf5730", "index": 8}, {"name": "file_data_store", "color": "#bf5730", "index": 8}, {"name": "qload", "color": "#bf5730", "index": 8}, {"name": "qsave", "color": "#bf5730", "index": 8}]}, {"name": "simdiag", "color": "#bf5730", "index": 8, "children": [{"name": "simdiag", "color": "#bf5730", "index": 8}]}, {"name": "version", "color": "#bf5730", "index": 8}]}
\ No newline at end of file
diff --git a/doc/apidoc/functions.rst b/doc/apidoc/functions.rst
index 8333a4127b..acacc97c80 100644
--- a/doc/apidoc/functions.rst
+++ b/doc/apidoc/functions.rst
@@ -162,8 +162,8 @@ Bloch-Redfield Master Equation
 Floquet States and Floquet-Markov Master Equation
 -------------------------------------------------
 
-.. automodule:: qutip.solve.floquet
-    :members: fmmesolve, 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_master_equation_steadystate, floquet_basis_transform, floquet_markov_mesolve
+.. automodule:: qutip.solver.floquet
+    :members: fmmesolve, fsesolve, FloquetBasis, FMESolver, floquet_tensor
 
 
 Stochastic Schrödinger Equation and Master Equation
diff --git a/doc/guide/dynamics/dynamics-floquet.rst b/doc/guide/dynamics/dynamics-floquet.rst
index c2affdd3c3..e3d0a3020c 100644
--- a/doc/guide/dynamics/dynamics-floquet.rst
+++ b/doc/guide/dynamics/dynamics-floquet.rst
@@ -93,7 +93,7 @@ Consider for example the case of a strongly driven two-level atom, described by
 In QuTiP we can define this Hamiltonian as follows:
 
 .. plot::
-   :context:
+   :context: close-figs
 
    >>> delta = 0.2 * 2*np.pi
    >>> eps0 = 1.0 * 2*np.pi
@@ -104,15 +104,17 @@ In QuTiP we can define this Hamiltonian as follows:
    >>> args = {'w': omega}
    >>> H = [H0, [H1, 'sin(w * t)']]
 
-The :math:`t=0` Floquet modes corresponding to the Hamiltonian :eq:`eq_driven_qubit` can then be calculated using the :func:`qutip.floquet.floquet_modes` function, which returns lists containing the Floquet modes and the quasienergies
+The :math:`t=0` Floquet modes corresponding to the Hamiltonian :eq:`eq_driven_qubit` can then be calculated using the :class:`qutip.FloquetBasis` class, which encapsulates the Floquet modes and the quasienergies:
 
 .. plot::
    :context:
 
    >>> T = 2*np.pi / omega
-   >>> f_modes_0, f_energies = floquet_modes(H, T, args)
+   >>> floquet_basis = FloquetBasis(H, T, args)
+   >>> f_energies = floquet_basis.e_quasi
    >>> f_energies # doctest: +NORMALIZE_WHITESPACE
    array([-2.83131212,  2.83131212])
+   >>> f_modes_0 = floquet_basis.mode(0)
    >>> f_modes_0 # doctest: +NORMALIZE_WHITESPACE
    [Quantum object: dims = [[2], [1]], shape = (2, 1), type = ket
    Qobj data =
@@ -123,26 +125,29 @@ The :math:`t=0` Floquet modes corresponding to the Hamiltonian :eq:`eq_driven_qu
    [[0.39993746+0.554682j]
     [0.72964231+0.j      ]]]
 
-For some problems interesting observations can be draw from the quasienergy levels alone. Consider for example the quasienergies for the driven two-level system introduced above as a function of the driving amplitude, calculated and plotted in the following example. For certain driving amplitudes the quasienergy levels cross. Since the quasienergies can be associated with the time-scale of the long-term dynamics due that the driving, degenerate quasienergies indicates a "freezing" of the dynamics (sometimes known as coherent destruction of tunneling).
+For some problems interesting observations can be draw from the quasienergy levels alone.
+Consider for example the quasienergies for the driven two-level system introduced above as a function of the driving amplitude, calculated and plotted in the following example.
+For certain driving amplitudes the quasienergy levels cross.
+Since the quasienergies can be associated with the time-scale of the long-term dynamics due that the driving, degenerate quasienergies indicates a "freezing" of the dynamics (sometimes known as coherent destruction of tunneling).
 
 .. plot::
    :context:
 
-   >>> delta = 0.2 * 2*np.pi
-   >>> eps0  = 0.0 * 2*np.pi
-   >>> omega = 1.0 * 2*np.pi
+   >>> delta = 0.2 * 2 * np.pi
+   >>> eps0  = 0.0 * 2 * np.pi
+   >>> omega = 1.0 * 2 * np.pi
    >>> A_vec = np.linspace(0, 10, 100) * omega
-   >>> T = (2*np.pi)/omega
+   >>> T = (2 * np.pi) / omega
    >>> tlist = np.linspace(0.0, 10 * T, 101)
-   >>> spsi0 = basis(2,0)
+   >>> spsi0 = basis(2, 0)
    >>> q_energies = np.zeros((len(A_vec), 2))
-   >>> H0 = delta/2.0 * sigmaz() - eps0/2.0 * sigmax()
+   >>> H0 = delta / 2.0 * sigmaz() - eps0 / 2.0 * sigmax()
    >>> args = {'w': omega}
    >>> for idx, A in enumerate(A_vec): # doctest: +SKIP
-   >>>   H1 = A/2.0 * sigmax() # doctest: +SKIP
-   >>>   H = [H0, [H1, lambda t, args: np.sin(args['w']*t)]] # doctest: +SKIP
-   >>>   f_modes, f_energies = floquet_modes(H, T, args, True) # doctest: +SKIP
-   >>>   q_energies[idx,:] = f_energies # doctest: +SKIP
+   >>>   H1 = A / 2.0 * sigmax() # doctest: +SKIP
+   >>>   H = [H0, [H1, lambda t, args: np.sin(args['w'] * t)]] # doctest: +SKIP
+   >>>   floquet_basis = FloquetBasis(H, T, args)
+   >>>   q_energies[idx,:] = floquet_basis.e_quasi # doctest: +SKIP
    >>> plt.figure() # doctest: +SKIP
    >>> plt.plot(A_vec/omega, q_energies[:,0] / delta, 'b', A_vec/omega, q_energies[:,1] / delta, 'r') # doctest: +SKIP
    >>> plt.xlabel(r'$A/\omega$') # doctest: +SKIP
@@ -150,12 +155,12 @@ For some problems interesting observations can be draw from the quasienergy leve
    >>> plt.title(r'Floquet quasienergies') # doctest: +SKIP
    >>> plt.show() # doctest: +SKIP
 
-Given the Floquet modes at :math:`t=0`, we obtain the Floquet mode at some later time :math:`t` using the function :func:`qutip.floquet.floquet_modes_t`:
+Given the Floquet modes at :math:`t=0`, we obtain the Floquet mode at some later time :math:`t` using :meth:`FloquetBasis.mode`:
 
 .. plot::
    :context: close-figs
 
-   >>> f_modes_t = floquet_modes_t(f_modes_0, f_energies, 2.5, H, T, args)
+   >>> f_modes_t = floquet_basis.mode(2.5)
    >>> f_modes_t # doctest: +SKIP
    [Quantum object: dims = [[2], [1]], shape = (2, 1), type = ket
    Qobj data =
@@ -166,24 +171,25 @@ Given the Floquet modes at :math:`t=0`, we obtain the Floquet mode at some later
    [[-0.37793106-0.00431336j]
     [-0.89630512+0.23191946j]]]
 
-The purpose of calculating the Floquet modes is to find the wavefunction solution to the original problem :eq:`eq_driven_qubit` given some initial state :math:`\left|\psi_0\right>`. To do that, we first need to decompose the initial state in the Floquet states, using the function :func:`qutip.floquet.floquet_state_decomposition`
+The purpose of calculating the Floquet modes is to find the wavefunction solution to the original problem :eq:`eq_driven_qubit` given some initial state :math:`\left|\psi_0\right>`.
+To do that, we first need to decompose the initial state in the Floquet states, using the function :meth:`FloquetBasis.to_floquet_basis`
 
 .. plot::
    :context:
 
    >>> psi0 = rand_ket(2)
-   >>> f_coeff = floquet_state_decomposition(f_modes_0, f_energies, psi0)
+   >>> f_coeff = floquet_basis.to_floquet_basis(psi0)
    >>> f_coeff # doctest: +SKIP
    [(-0.645265993068382+0.7304552549315746j),
    (0.15517002114250228-0.1612116102238258j)]
 
-and given this decomposition of the initial state in the Floquet states we can easily evaluate the wavefunction that is the solution to :eq:`eq_driven_qubit` at an arbitrary time :math:`t` using the function :func:`qutip.floquet.floquet_wavefunction_t`
+and given this decomposition of the initial state in the Floquet states we can easily evaluate the wavefunction that is the solution to :eq:`eq_driven_qubit` at an arbitrary time :math:`t` using the function :meth:`FloquetBasis.from_floquet_basis`:
 
 .. plot::
    :context:
 
    >>> t = 10 * np.random.rand()
-   >>> psi_t = floquet_wavefunction_t(f_modes_0, f_energies, f_coeff, t, H, T, args)
+   >>> psi_t = floquet_basis.from_floquet_basis(f_coeff, t)
 
 The following example illustrates how to use the functions introduced above to calculate and plot the time-evolution of :eq:`eq_driven_qubit`.
 
@@ -194,7 +200,9 @@ The following example illustrates how to use the functions introduced above to c
 Pre-computing the Floquet modes for one period
 ----------------------------------------------
 
-When evaluating the Floquet states or the wavefunction at many points in time it is useful to pre-compute the Floquet modes for the first period of the driving with the required resolution. In QuTiP the function :func:`qutip.floquet.floquet_modes_table` calculates a table of Floquet modes which later can be used together with the function :func:`qutip.floquet.floquet_modes_t_lookup` to efficiently lookup the Floquet mode at an arbitrary time. The following example illustrates how the example from the previous section can be solved more efficiently using these functions for pre-computing the Floquet modes.
+When evaluating the Floquet states or the wavefunction at many points in time it is useful to pre-compute the Floquet modes for the first period of the driving with the required times.
+The list of times to pre-compute modes for may be passed to :class:`FloquetBasis` using `precompute=tlist`, and then `:meth:`FloquetBasis.from_floquet_basis` and :meth:`FloquetBasis.to_floquet_basis` can be used to efficiently retrieve the wave function at the pre-computed times.
+The following example illustrates how the example from the previous section can be solved more efficiently using these functions for pre-computing the Floquet modes:
 
 .. plot:: guide/scripts/floquet_ex2.py
    :width: 4.0in
@@ -231,9 +239,9 @@ For any coupling operator :math:`q` (given by the user) the matrix elements in t
 From the matrix elements and the spectral density :math:`J(\omega)`, the decay rate :math:`\gamma_{\alpha \beta k}` is defined:
 
 .. math::
-    \gamma_{\alpha \beta k} = 2 \pi \Theta(\Delta_{\alpha \beta k}) J(\Delta_{\alpha \beta k}) | X_{\alpha \beta k}|^2
+    \gamma_{\alpha \beta k} = 2 \pi J(\Delta_{\alpha \beta k}) | X_{\alpha \beta k}|^2
 
-where :math:`\Theta` is the Heaviside function. The master equation is further simplified by the RWA, which makes the following matrix useful:
+The master equation is further simplified by the RWA, which makes the following matrix useful:
 
 .. math::
     A_{\alpha \beta} = \sum_{k = -\infty}^\infty [\gamma_{\alpha \beta k} + n_{th}(|\Delta_{\alpha \beta k}|)(\gamma_{\alpha \beta k} + \gamma_{\alpha \beta -k})
@@ -263,7 +271,7 @@ The noise spectral-density function of the environment is implemented as a Pytho
 
     gamma1 = 0.1
     def noise_spectrum(omega):
-        return 0.5 * gamma1 * omega/(2*pi)
+        return (omega>0) * 0.5 * gamma1 * omega/(2*pi)
 
 The other parameters are similar to the :func:`qutip.mesolve` and :func:`qutip.mcsolve`, and the same format for the return value is used :class:`qutip.solve.solver.Result`. The following example extends the example studied above, and uses :func:`qutip.floquet.fmmesolve` to introduce dissipation into the calculation
 
@@ -271,7 +279,7 @@ The other parameters are similar to the :func:`qutip.mesolve` and :func:`qutip.m
    :width: 4.0in
    :include-source:
 
-Alternatively, we can let the :func:`qutip.floquet.fmmesolve` function transform the density matrix at each time step back to the computational basis, and calculating the expectation values for us, by using:
+Finally, :func:`qutip.solver.floquet.fmmesolve`  always expects the ``e_ops`` to be specified in the laboratory basis (as for other solvers) and we can calculate expectation values using:
 
-    output = fmmesolve(H, psi0, tlist, [sigmax()], [num(2)], [noise_spectrum], T, args, floquet_basis=False)
+    output = fmmesolve(H, psi0, tlist, [sigmax()], e_ops=[num(2)], spectra_cb=[noise_spectrum], T=T, args=args)
     p_ex = output.expect[0]
diff --git a/doc/guide/scripts/floquet_ex0.py b/doc/guide/scripts/floquet_ex0.py
deleted file mode 100644
index 0b3165b4f2..0000000000
--- a/doc/guide/scripts/floquet_ex0.py
+++ /dev/null
@@ -1,31 +0,0 @@
-import numpy as np
-from matplotlib import pyplot
-import qutip
-
-delta = 0.2 * 2*np.pi
-eps0  = 0.0 * 2*np.pi
-omega = 1.0 * 2*np.pi
-A_vec = np.linspace(0, 10, 100) * omega
-T      = 2*np.pi/omega
-tlist  = np.linspace(0.0, 10 * T, 101)
-psi0   = qutip.basis(2, 0)
-
-q_energies = np.zeros((len(A_vec), 2))
-
-H0 = delta/2.0 * qutip.sigmaz() - eps0/2.0 * qutip.sigmax()
-args = {'w': omega}
-for idx, A in enumerate(A_vec):
-    H1 = A/2.0 * qutip.sigmax()
-    H = [H0, [H1, lambda t, w: np.sin(w*t)]]
-    f_modes,f_energies = qutip.floquet_modes(H, T, args, True)
-    q_energies[idx,:] = f_energies
-
-# plot the results
-pyplot.plot(
-    A_vec/omega, np.real(q_energies[:, 0]) / delta, 'b',
-    A_vec/omega, np.real(q_energies[:, 1]) / delta, 'r',
-)
-pyplot.xlabel(r'$A/\omega$')
-pyplot.ylabel(r'Quasienergy / $\Delta$')
-pyplot.title(r'Floquet quasienergies')
-pyplot.show()
diff --git a/doc/guide/scripts/floquet_ex1.py b/doc/guide/scripts/floquet_ex1.py
index 7afea3efa2..3ba5646c17 100644
--- a/doc/guide/scripts/floquet_ex1.py
+++ b/doc/guide/scripts/floquet_ex1.py
@@ -14,18 +14,18 @@
 H0 = - delta/2.0 * qutip.sigmax() - eps0/2.0 * qutip.sigmaz()
 H1 = A/2.0 * qutip.sigmaz()
 args = {'w': omega}
-H = [H0, [H1, lambda t,args: np.sin(args['w'] * t)]]
+H = [H0, [H1, lambda t, w: np.sin(w * t)]]
 
-# find the floquet modes for the time-dependent hamiltonian
-f_modes_0,f_energies = qutip.floquet_modes(H, T, args)
+# Create the floquet system for the time-dependent hamiltonian
+floquetbasis = qutip.FloquetBasis(H, T, args)
 
 # decompose the inital state in the floquet modes
-f_coeff = qutip.floquet_state_decomposition(f_modes_0, f_energies, psi0)
+f_coeff = floquetbasis.to_floquet_basis(psi0)
 
-# calculate the wavefunctions using the from the floquet modes
+# calculate the wavefunctions using the from the floquet modes coefficients
 p_ex = np.zeros(len(tlist))
 for n, t in enumerate(tlist):
-    psi_t = qutip.floquet_wavefunction_t(f_modes_0, f_energies, f_coeff, t, H, T, args)
+    psi_t = floquetbasis.from_floquet_basis(f_coeff, t)
     p_ex[n] = qutip.expect(qutip.num(2), psi_t)
 
 # For reference: calculate the same thing with mesolve
diff --git a/doc/guide/scripts/floquet_ex2.py b/doc/guide/scripts/floquet_ex2.py
index fd652bc76b..823dec0b7d 100644
--- a/doc/guide/scripts/floquet_ex2.py
+++ b/doc/guide/scripts/floquet_ex2.py
@@ -13,20 +13,18 @@
 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, args: np.sin(args['w'] * t)]]
+H = [H0, [H1, lambda t, w: np.sin(w * t)]]
 
-# find the floquet modes for the time-dependent hamiltonian        
-f_modes_0,f_energies = qutip.floquet_modes(H, T, args)
+# find the floquet modes for the time-dependent hamiltonian
+floquetbasis = qutip.FloquetBasis(H, T, args, precompute=tlist)
 
 # decompose the inital state in the floquet modes
-f_coeff = qutip.floquet_state_decomposition(f_modes_0, f_energies, psi0)
+f_coeff = floquetbasis.to_floquet_basis(psi0)
 
-# calculate the wavefunctions using the from the floquet modes
-f_modes_table_t = qutip.floquet_modes_table(f_modes_0, f_energies, tlist, H, T, args)
+# calculate the wavefunctions using the from the floquet modes coefficients
 p_ex = np.zeros(len(tlist))
 for n, t in enumerate(tlist):
-    f_modes_t = qutip.floquet_modes_t_lookup(f_modes_table_t, t, T)
-    psi_t     = qutip.floquet_wavefunction(f_modes_t, f_energies, f_coeff, t)
+    psi_t = floquetbasis.from_floquet_basis(f_coeff, t)
     p_ex[n] = qutip.expect(qutip.num(2), psi_t)
 
 # For reference: calculate the same thing with mesolve
diff --git a/doc/guide/scripts/floquet_ex3.py b/doc/guide/scripts/floquet_ex3.py
index 06f3fedb1d..d6f2b75b63 100644
--- a/doc/guide/scripts/floquet_ex3.py
+++ b/doc/guide/scripts/floquet_ex3.py
@@ -7,36 +7,34 @@
 A     = 0.25 * 2*np.pi
 omega = 1.0 * 2*np.pi
 T      = 2*np.pi / omega
-tlist  = np.linspace(0.0, 20 * T, 101)
+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,args: np.sin(args['w'] * t)]]
+H = [H0, [H1, lambda t, w: np.sin(w * t)]]
 
 # noise power spectrum
 gamma1 = 0.1
 def noise_spectrum(omega):
-    return 0.5 * gamma1 * omega/(2*np.pi)
-
-# find the floquet modes for the time-dependent hamiltonian        
-f_modes_0, f_energies = qutip.floquet_modes(H, T, args)
-
-# precalculate mode table
-f_modes_table_t = qutip.floquet_modes_table(
-    f_modes_0, f_energies, np.linspace(0, T, 500 + 1), H, T, args,
-)
+    return (omega>0) * 0.5 * gamma1 * omega/(2*np.pi)
 
 # solve the floquet-markov master equation
-output = qutip.fmmesolve(H, psi0, tlist, [qutip.sigmax()], [], [noise_spectrum], T, args)
+output = qutip.fmmesolve(
+    H, psi0, tlist, [qutip.sigmax()],
+    spectra_cb=[noise_spectrum], T=T,
+    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_modes_t = qutip.floquet_modes_t_lookup(f_modes_table_t, t, T)
-    f_states_t = qutip.floquet_states(f_modes_t, f_energies, t)
-    p_ex[idx] = qutip.expect(qutip.num(2), output.states[idx].transform(f_states_t, True))
+    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,
diff --git a/qutip/solve/__init__.py b/qutip/solve/__init__.py
index 2e7cf99238..cd275059f4 100644
--- a/qutip/solve/__init__.py
+++ b/qutip/solve/__init__.py
@@ -1,4 +1,3 @@
-from .floquet import *
 from .mcsolve import *
 from .mesolve import *
 from . import nonmarkov
diff --git a/qutip/solve/floquet.py b/qutip/solve/floquet.py
deleted file mode 100644
index 0eec45e813..0000000000
--- a/qutip/solve/floquet.py
+++ /dev/null
@@ -1,1058 +0,0 @@
-__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
-import warnings
-from copy import copy
-from numpy import angle, pi, exp, sqrt
-from types import FunctionType
-from .. import (
-    Qobj, unstacked_index, stack_columns, unstack_columns, projection, expect,
-)
-from ..core import data as _data
-from .sesolve import sesolve
-from ._steadystate_v4 import steadystate
-from .solver import SolverOptions
-from ..solver.propagator import propagator
-from .solver import Result, _solver_safety_check
-from ..utilities import n_thermal
-
-
-
-def floquet_modes(H, T, args=None, sort=False, U=None, options=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`.
-
-    options : :class:`qutip.solver.Options`
-        options for the ODE solver. For the propagator U.
-
-    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=args, options=copy(options))
-
-    # 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.array(evecs[:, o]).T,
-                       dims=new_dims) for o in order]
-
-    return kets_order, e_quasi[order]
-
-
-def floquet_modes_t(f_modes_0, f_energies, t, H, T, args=None,
-                    options=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.
-
-    options : :class:`qutip.solver.Options`
-        options for the ODE solver. For the propagator.
-
-    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, options=copy(options))
-
-        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
-
-
-def floquet_modes_table(f_modes_0, f_energies, tlist, H, T, args=None,
-                        options=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
-
-    options : :class:`qutip.solver.Options`
-        options for the ODE solver.
-
-    Returns
-    -------
-
-    output : nested list
-
-        A nested list of Floquet modes as kets for each time in `tlist`
-
-    """
-    options = copy(options) or SolverOptions()
-    # truncate tlist to the driving period
-    tlist_period = tlist[np.where(tlist <= T)]
-
-    f_modes_table_t = [[] for t in tlist_period]
-
-    opt = SolverOptions()
-
-    for n, f_mode in enumerate(f_modes_0):
-        output = sesolve(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
-
-
-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))]
-
-
-def floquet_states_t(f_modes_0, f_energies, t, H, T, args=None,
-                     options=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.
-
-    options : :class:`qutip.solver.Options`
-        options for the ODE solver.
-
-    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,
-                                options=options)
-    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(1, len(f_energies))],
-        start=f_modes_t[0] * exp(-1j * f_energies[0] * t) * f_coeff[0]
-    )
-
-
-def floquet_wavefunction_t(f_modes_0, f_energies, f_coeff, t, H, T, args=None,
-                           options=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,
-                                  options=options)
-    return sum(
-        [f_states_t[i] * f_coeff[i] for i in np.arange(1, len(f_energies))],
-        start = f_states_t[0] * f_coeff[0]
-    )
-
-
-def floquet_state_decomposition(f_states, f_energies, psi):
-    r"""
-    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 [state.dag() * psi for state in f_states]
-
-
-
-def fsesolve(H, psi0, tlist, e_ops=[], T=None, args={}, Tsteps=100,
-             options_modes=None):
-    """
-    Solve the Schrodinger equation using the Floquet formalism.
-
-    Parameters
-    ----------
-
-    H : :class:`qutip.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.
-
-    options_modes : :class:`qutip.solver.Options`
-        options for the ODE solver.
-
-    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]
-
-    if options_modes is None:
-        options_modes_table = SolverOptions()
-    else:
-        options_modes_table = options_modes
-
-    # find the floquet modes for the time-dependent hamiltonian
-    f_modes_0, f_energies = floquet_modes(H, T, args,
-                                          options=options_modes)
-
-    # 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,
-                                          options=options_modes_table)
-
-    # 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.
-
-    .. note :
-        The number of integration steps (for calculating X) within one period
-        is set to 20 * kmax.
-
-    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.
-
-    kmax : 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).
-
-    options : :class:`qutip.solver.Options`
-        options for the ODE solver.
-
-    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))
-
-    # time steps for integration of coupling operator
-    nT = int(np.max([20 * kmax, 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 = np.hstack([f.full() for f in floquet_modes_t_lookup(
-            f_modes_table_t, t, T)])
-        FF = f_modes_t.T.conj() @ c_op.full() @ f_modes_t
-        phi = exp(-1j * np.arange(-kmax, kmax+1) * omega * t)
-        X += (dT / T) * np.einsum("ij,k->ijk", FF, phi)
-
-    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])
-
-    Rdata_lil = scipy.sparse.lil_matrix((N * N, N * N), dtype=complex)
-
-    AsumList = [np.sum(A, axis=1) for A in Alist]
-
-    for k in range(len(Alist)):
-        for i in range(N):
-            Rdata_lil[i+N*i, i+N*i] -= -Alist[k][i, i] + AsumList[k][i]
-            for j in range(i+1, N):
-                Rdata_lil[i+N*i, j+N*j] += Alist[k][j, i]
-                Rdata_lil[j+N*j, i+N*i] += Alist[k][i, j]
-                a_term = -(1/2)*(AsumList[k][i] + AsumList[k][j])
-                Rdata_lil[i+N*j, i+N*j] += a_term
-                Rdata_lil[j+N*i, j+N*i] += a_term
-
-    return Qobj(Rdata_lil, dims=[[N, N], [N, N]])
-
-
-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, rho0, tlist, e_ops, options=None, floquet_basis=True,
-    f_modes_0=None, f_modes_table_t=None, f_energies=None, T=None,
-):
-    """
-    Solve the dynamics for the system using the Floquet-Markov master equation.
-
-    .. note::
-
-        It is important to understand in which frame and basis the results
-        are returned here.
-
-    Parameters
-    ----------
-
-    R : array
-        The Floquet-Markov master equation tensor `R`.
-
-    rho0 : :class:`qutip.qobj`
-        Initial density matrix.  If ``f_modes_0`` is not passed, this density
-        matrix is assumed to be in the Floquet picture.
-
-    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.
-
-    options : :class:`qutip.solver.Options`
-        options for the ODE solver.
-
-    floquet_basis: bool, True
-        If ``True``, states and expectation values will be returned in the
-        Floquet basis.  If ``False``, a transformation will be made to the
-        computational basis; this will be in the lab frame if
-        ``f_modes_table``, ``T` and ``f_energies`` are all supplied, or the
-        interaction picture (defined purely be f_modes_0) if they are not.
-
-    f_modes_0 : list of :class:`qutip.qobj` (kets), optional
-        A list of initial Floquet modes, used to transform the given starting
-        density matrix into the Floquet basis.  If this is not passed, it is
-        assumed that ``rho`` is already in the Floquet basis.
-
-    f_modes_table_t : nested list of :class:`qutip.qobj` (kets), optional
-        A lookup-table of Floquet modes at times precalculated by
-        :func:`qutip.floquet.floquet_modes_table`.  Necessary if
-        ``floquet_basis`` is ``False`` and the transformation should be made
-        back to the lab frame.
-
-    f_energies : array_like of float, optional
-        The precalculated Floquet quasienergies.  Necessary if
-        ``floquet_basis`` is ``False`` and the transformation should be made
-        back to the lab frame.
-
-    T : float, optional
-        The time period of driving.  Necessary if ``floquet_basis`` is
-        ``False`` and the transformation should be made back to the lab frame.
-
-    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`.
-    """
-    opt = options or SolverOptions()
-    rho0 = rho0.proj() if rho0.isket else rho0
-
-    # Prepare output object.
-    dt = tlist[1] - tlist[0]
-    output = Result()
-    output.solver = "fmmesolve"
-    output.times = tlist
-    if isinstance(e_ops, FunctionType):
-        expt_callback = True
-        store_states = opt['store_states'] or False
-    else:
-        expt_callback = False
-        try:
-            e_ops = list(e_ops)
-        except TypeError:
-            raise TypeError("`e_ops` must be iterable or a function") from None
-        n_expt_op = len(e_ops)
-        if n_expt_op == 0:
-            store_states = True
-        else:
-            output.expect = []
-            output.num_expect = n_expt_op
-            for op in e_ops:
-                dtype = np.float64 if op.isherm else np.complex128
-                output.expect.append(np.zeros(len(tlist), dtype=dtype))
-        store_states = opt['store_states'] or (n_expt_op == 0)
-    if store_states:
-        output.states = []
-
-    # Choose which frame transformations should be done on the initial and
-    # evolved states.
-    lab_lookup = [f_modes_table_t, f_energies, T]
-    if (
-        any(x is None for x in lab_lookup)
-        and not all(x is None for x in lab_lookup)
-    ):
-        warnings.warn(
-            "if transformation back to the computational basis in the lab"
-            "frame is desired, all of `f_modes_t`, `f_energies` and `T` must"
-            "be supplied."
-        )
-        f_modes_table_t = f_energies = T = None
-
-    # Initial state.
-    if f_modes_0 is not None:
-        rho0 = rho0.transform(f_modes_0)
-
-    # Evolved states.
-    if floquet_basis:
-        def transform(rho, t):
-            return rho
-    elif f_modes_table_t is not None:
-        # Lab frame, computational basis.
-        def transform(rho, t):
-            f_modes_t = floquet_modes_t_lookup(f_modes_table_t, t, T)
-            f_states_t = floquet_states(f_modes_t, f_energies, t)
-            return rho.transform(f_states_t, True)
-    elif f_modes_0 is not None:
-        # Interaction picture, computational basis.
-        def transform(rho, t):
-            return rho.transform(f_modes_0, False)
-    else:
-        raise ValueError(
-            "cannot transform out of the Floquet basis without some knowledge "
-            "of the Floquet modes.  Pass `f_modes_0`, or all of `f_modes_t`, "
-            "`f_energies` and `T`."
-        )
-
-    #
-    # setup integrator
-    #
-    initial_vector = stack_columns(rho0.full())
-    r = scipy.integrate.ode(_wrap_matmul)
-    r.set_f_params(R.data)
-    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])
-
-    # Main evolution loop.
-    for t_idx, t in enumerate(tlist):
-        if not r.successful():
-            break
-
-        rho = transform(Qobj(unstack_columns(r.y), rho0.dims), t)
-
-        if expt_callback:
-            e_ops(t, rho)
-        else:
-            for m, e_op in enumerate(e_ops):
-                output.expect[m][t_idx] = expect(e_op, rho)
-        if store_states:
-            output.states.append(rho)
-        r.integrate(r.t + dt)
-    return output
-
-
-def _wrap_matmul(t, state, operator):
-    return _data.matmul(operator, _data.dense.fast_from_numpy(state),
-                        dtype=_data.Dense).as_ndarray()
-
-# -----------------------------------------------------------------------------
-# Solve the Floquet-Markov master equation
-#
-#
-
-
-def fmmesolve(H, rho0, tlist, c_ops=[], e_ops=[], spectra_cb=[], T=None,
-              args={}, options=SolverOptions(), floquet_basis=True, kmax=5,
-              _safe_mode=True, options_modes=None):
-    """
-    Solve the dynamics for the system using the Floquet-Markov master equation.
-
-    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 # doctest: +SKIP
-        >>> h = 6.626e-34 # doctest: +SKIP
-        >>> kB = 1.38e-23 # doctest: +SKIP
-        >>> args['w_th'] = temperature * (kB / h) * 2 * pi * 1e-9 \
-            #doctest: +SKIP
-
-    options : :class:`qutip.solver.Options`
-        options for the ODE solver. For solving the master equation.
-
-    floquet_basis : bool
-        Will return results in Floquet basis or computational basis
-        (optional).
-
-    k_max : int
-        The truncation of the number of sidebands (default 5).
-
-    options_modes : :class:`qutip.solver.Options`
-        options for the ODE solver. For computing Floquet modes.
-
-    Returns
-    -------
-
-    output : :class:`qutip.solver.Result`
-
-        An instance of the class :class:`qutip.solver.Result`, which contains
-        either an *array* of expectation values for the times specified
-        by `tlist`.
-    """
-
-    if _safe_mode:
-        _solver_safety_check(H, rho0, c_ops, e_ops, args)
-
-    if options_modes is None:
-        options_modes_table = SolverOptions()
-    else:
-        options_modes_table = options_modes
-
-    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)
-
-    if len(spectra_cb) != len(c_ops):
-        raise ValueError("Length of c_ops and spectra_cb don't match.")
-
-    f_modes_0, f_energies = floquet_modes(H, T, args,
-                                          options=options_modes)
-
-    f_modes_table_t = floquet_modes_table(f_modes_0, f_energies,
-                                          np.linspace(0, T, 500 + 1),
-                                          H, T, args,
-                                          options=options_modes_table)
-
-    # get w_th from args if it exists
-    if 'w_th' in args:
-        w_th = args['w_th']
-    else:
-        w_th = 0
-
-    # floquet-markov master equation tensor
-    R = 0
-    # loop over input c_ops and spectra_cb, calculate one R for each set
-    for c_op, spectrum in zip(c_ops, spectra_cb):
-        # calculate the rate-matrices for the floquet-markov master equation
-        Delta, X, Gamma, Amat = floquet_master_equation_rates(
-            f_modes_0, f_energies, c_op, H, T, args, spectrum,
-            w_th, kmax, f_modes_table_t)
-
-        # calculate temporary floquet-markov master equation tensor
-        R += floquet_master_equation_tensor(Amat, f_energies)
-
-    return floquet_markov_mesolve(R, rho0, tlist, e_ops,
-                                  options=options,
-                                  floquet_basis=floquet_basis,
-                                  f_modes_0=f_modes_0,
-                                  f_modes_table_t=f_modes_table_t,
-                                  T=T,
-                                  f_energies=f_energies)
diff --git a/qutip/solver/__init__.py b/qutip/solver/__init__.py
index 226e76021a..1a4af81b70 100644
--- a/qutip/solver/__init__.py
+++ b/qutip/solver/__init__.py
@@ -6,4 +6,6 @@
 from .scattering import *
 from .correlation import *
 from .spectrum import *
+from .floquet import *
+from .floquet_bwcomp import *
 from .steadystate import *
diff --git a/qutip/solver/floquet.py b/qutip/solver/floquet.py
new file mode 100644
index 0000000000..25fca8aff7
--- /dev/null
+++ b/qutip/solver/floquet.py
@@ -0,0 +1,932 @@
+__all__ = [
+    "FloquetBasis",
+    "floquet_tensor",
+    "fsesolve",
+    "fmmesolve",
+    "FMESolver",
+]
+
+import numpy as np
+from qutip.core import data as _data
+from qutip import Qobj, QobjEvo
+from .propagator import Propagator
+from .mesolve import MESolver
+from .solver_base import Solver
+from .integrator import Integrator
+from .result import Result
+from time import time
+from ..ui.progressbar import progess_bars
+
+
+class FloquetBasis:
+    """
+    Utility to compute floquet modes and states.
+
+    Attributes
+    ----------
+    U : :class:`Propagator`
+        The propagator of the Hamiltonian over one period.
+
+    evecs : :class:`qutip.data.Data`
+        Matrix where each column is an initial Floquet mode.
+
+    e_quasi : np.ndarray[float]
+        The quasi energies of the Hamiltonian.
+    """
+
+    def __init__(
+        self,
+        H,
+        T,
+        args=None,
+        options=None,
+        sparse=False,
+        sort=True,
+        precompute=None,
+    ):
+        """
+        Parameters
+        ----------
+        H : :class:`Qobj`, :class:`QobjEvo`, QobjEvo compatible format.
+            System Hamiltonian, with period `T`.
+
+        T : float
+            Period of the Hamiltonian.
+
+        args : None / *dictionary*
+            dictionary of parameters for time-dependent Hamiltonians and
+            collapse operators.
+
+        options : dict [None]
+            Options used by sesolve to compute the floquet modes.
+
+        sparse : bool [False]
+            Whether to use the sparse eigen solver when computing the
+            quasi-energies.
+
+        sort : bool [True]
+            Whether to sort the quasi-energies.
+
+        precompute : list [None]
+            If provided, a list of time at which to store the propagators
+            for later use when computing modes and states. Default is
+            ``linspace(0, T, 101)`` corresponding to the default integration
+            steps used for the floquet tensor computation.
+        """
+        if not T > 0:
+            raise ValueError("The period need to be a positive number.")
+        self.T = T
+        if precompute is not None:
+            tlist = np.unique(np.atleast_1d(precompute) % self.T)
+            memoize = len(tlist)
+            if tlist[0] != 0:
+                memoize += 1
+            if tlist[-1] != T:
+                memoize += 1
+        else:
+            # Default computation
+            tlist = np.linspace(0, T, 101)
+            memoize = 101
+        self.U = Propagator(H, args=args, options=options, memoize=memoize)
+        for t in tlist:
+            # Do the evolution by steps to save the intermediate results.
+            self.U(t)
+        U_T = self.U(self.T)
+        if not sparse and isinstance(U_T.data, _data.CSR):
+            U_T = U_T.to("Dense")
+        evals, evecs = _data.eigs(U_T.data)
+        e_quasi = -np.angle(evals) / T
+        if sort:
+            perm = np.argsort(e_quasi)
+            self.evecs = _data.permute.indices(evecs, col_perm=np.argsort(perm))
+            self.e_quasi = e_quasi[perm]
+        else:
+            self.evecs = evecs
+            self.e_quasi = e_quasi
+
+    def _as_ketlist(self, kets_mat):
+        """
+        Split the Data array in a list of kets.
+        """
+        dims = [self.U(0).dims[0], [1]]
+        return [
+            Qobj(ket, dims=dims, type="ket")
+            for ket in _data.split_columns(kets_mat)
+        ]
+
+    def mode(self, t, data=False):
+        """
+        Calculate the Floquet modes at time ``t``.
+
+        Parameters
+        ----------
+        t : float
+            The time for which to evaluate the Floquet mode.
+
+        data : bool [False]
+            Whether to return the states as a single data matrix or a list of
+            ket states.
+
+        Returns
+        -------
+        output : list[:class:`Qobj`], :class:`qutip.data.Data`
+            A list of Floquet states for the time ``t`` or the states as column
+            in a single matrix.
+        """
+        t = t % self.T
+        if t == 0.0:
+            kets_mat = self.evecs
+        else:
+            U = self.U(t).data
+            phases = _data.diag(np.exp(1j * t * self.e_quasi))
+            kets_mat = U @ self.evecs @ phases
+        if data:
+            return kets_mat
+        else:
+            return self._as_ketlist(kets_mat)
+
+    def state(self, t, data=False):
+        """
+        Evaluate the floquet states at time t.
+
+        Parameters
+        ----------
+        t : float
+            The time for which to evaluate the Floquet states.
+
+        data : bool [False]
+            Whether to return the states as a single data matrix or a list of
+            ket states.
+
+        Returns
+        -------
+        output : list[:class:`Qobj`], :class:`qutip.data.Data`
+            A list of Floquet states for the time ``t`` or the states as column
+            in a single matrix.
+        """
+        if t:
+            phases = _data.diag(np.exp(-1j * t * self.e_quasi))
+            states_mat = self.mode(t, True) @ phases
+        else:
+            states_mat = self.evecs
+        if data:
+            return states_mat
+        else:
+            return self._as_ketlist(states_mat)
+
+    def from_floquet_basis(self, floquet_basis, t=0):
+        """
+        Transform a ket or density matrix from the Floquet basis at time ``t``
+        to the lab basis.
+
+        Parameters
+        ----------
+        floquet_basis : :class:`Qobj`, :class:`qutip.data.Data`
+            Initial state in the Floquet basis at time ``t``. May be either a
+            ket or density matrix.
+
+        t : float [0]
+            The time at which to evaluate the Floquet states.
+
+        Returns
+        -------
+        output : :class:`Qobj`, :class:`qutip.data.Data`
+            The state in the lab basis. The return type is the same as the type
+            of the input state.
+        """
+        is_Qobj = isinstance(floquet_basis, Qobj)
+        if is_Qobj:
+            dims = floquet_basis.dims
+            floquet_basis = floquet_basis.data
+            if dims[0] != self.U(0).dims[1]:
+                raise ValueError(
+                    "Dimensions of the state do not match the Hamiltonian"
+                )
+
+        state_mat = self.state(t, True)
+        lab_basis = state_mat @ floquet_basis
+        if floquet_basis.shape[1] != 1:
+            lab_basis = lab_basis @ state_mat.adjoint()
+
+        if is_Qobj:
+            return Qobj(lab_basis, dims=dims)
+        return lab_basis
+
+    def to_floquet_basis(self, lab_basis, t=0):
+        """
+        Transform a ket or density matrix in the lab basis
+        to the Floquet basis at time ``t``.
+
+        Parameters
+        ----------
+        lab_basis : :class:`Qobj`, :class:`qutip.data.Data`
+            Initial state in the lab basis.
+
+        t : float [0]
+            The time at which to evaluate the Floquet states.
+
+        Returns
+        -------
+        output : :class:`Qobj`, :class:`qutip.data.Data`
+            The state in the Floquet basis. The return type is the same as the
+            type of the input state.
+        """
+        is_Qobj = isinstance(lab_basis, Qobj)
+        if is_Qobj:
+            dims = lab_basis.dims
+            lab_basis = lab_basis.data
+            if dims[0] != self.U(0).dims[1]:
+                raise ValueError(
+                    "Dimensions of the state do not match the Hamiltonian"
+                )
+
+        state_mat = self.state(t, True)
+        floquet_basis = state_mat.adjoint() @ lab_basis
+        if lab_basis.shape[1] != 1:
+            floquet_basis = floquet_basis @ state_mat
+
+        if is_Qobj:
+            return Qobj(floquet_basis, dims=dims)
+        return floquet_basis
+
+
+def _floquet_delta_tensor(f_energies, kmax, T):
+    """
+    Floquet-Markov master equation X matrices.
+
+    Parameters
+    ----------
+    f_energies : np.ndarray
+        The Floquet energies.
+
+    kmax : int
+        The truncation of the number of sidebands (default 5).
+
+    T : float
+        The period of the time-dependence of the Hamiltonian.
+
+    Returns
+    -------
+    delta : np.ndarray
+        Floquet delta tensor.
+    """
+    delta = np.subtract.outer(f_energies, f_energies)
+    return np.add.outer(delta, np.arange(-kmax, kmax + 1) * (2 * np.pi / T))
+
+
+def _floquet_X_matrices(floquet_basis, c_ops, kmax, ntimes=100):
+    """
+    Floquet-Markov master equation X matrices.
+
+    Parameters
+    ----------
+    floquet_basis : :class:`FloquetBasis`
+        The system Hamiltonian wrapped in a FloquetBasis object.
+
+    c_ops : list of :class:`Qobj`
+        The collapse operators describing the dissipation.
+
+    kmax : int
+        The truncation of the number of sidebands (default 5).
+
+    ntimes : int [100]
+        The number of integration steps (for calculating X) within one period.
+
+    Returns
+    -------
+    X : list of dict of :class:`qutip.data.Data`
+        A dict of the sidebands ``k`` for the X matrices of each c_ops
+    """
+    T = floquet_basis.T
+    N = floquet_basis.U(0).shape[0]
+    omega = (2 * np.pi) / T
+    tlist = np.linspace(T / ntimes, T, ntimes)
+    ks = np.arange(-kmax, kmax + 1)
+    out = {k: [_data.csr.zeros(N, N)] * len(c_ops) for k in ks}
+
+    for t in tlist:
+        mode = floquet_basis.mode(t, data=True)
+        FFs = [mode.adjoint() @ c_op.data @ mode for c_op in c_ops]
+        for k, phi in zip(ks, np.exp(-1j * ks * omega * t) / ntimes):
+            out[k] = [
+                _data.add(prev, new, phi) for prev, new in zip(out[k], FFs)
+            ]
+
+    return [{k: out[k][i] for k in ks} for i in range(len(c_ops))]
+
+
+def _floquet_gamma_matrices(X, delta, J_cb):
+    """
+    Floquet-Markov master equation gamma matrices.
+
+    Parameters
+    ----------
+    X : list of dict of :class:`qutip.data.Data`
+        Floquet X matrices created by :func:`_floquet_X_matrices`.
+
+    delta : np.ndarray
+        Floquet delta tensor created by :func:`_floquet_delta_tensor`.
+
+    J_cb : list of callables
+        A list callback functions that compute the noise power spectrum as
+        a function of frequency. The list should contain one callable for each
+        collapse operator `c_op`, in the same order as the elements of `X`.
+        Each callable should accept a numpy array of frequencies and return a
+        numpy array of corresponding noise power.
+
+    Returns
+    -------
+    gammas : dict of :class:`qutip.data.Data`
+        A dict mapping the sidebands ``k`` to their gamma matrices.
+    """
+    N = delta.shape[0]
+    kmax = (delta.shape[2] - 1) // 2
+    gamma = {k: _data.csr.zeros(N, N) for k in range(-kmax, kmax + 1, 1)}
+
+    for X_c_op, sp in zip(X, J_cb):
+        response = sp(delta) * ((2 + 0j) * np.pi)
+        response = [
+            _data.Dense(response[:, :, k], copy=False)
+            for k in range(2 * kmax + 1)
+        ]
+        for k in range(-kmax, kmax + 1, 1):
+            gamma[k] = _data.add(
+                gamma[k],
+                _data.multiply(
+                    _data.multiply(X_c_op[k].conj(), X_c_op[k]),
+                    response[k + kmax],
+                ),
+            )
+    return gamma
+
+
+def _floquet_A_matrix(delta, gamma, w_th):
+    """
+    Floquet-Markov master equation rate matrix.
+
+    Parameters
+    ----------
+    delta : np.ndarray
+        Floquet delta tensor created by :func:`_floquet_delta_tensor`.
+
+    gamma : dict of :class:`qutip.data.Data`
+        Floquet gamma matrices created by :func:`_floquet_gamma_matrices`.
+
+    w_th : float
+        The temperature in units of frequency.
+    """
+    kmax = (delta.shape[2] - 1) // 2
+
+    if w_th > 0.0:
+        deltap = np.copy(delta)
+        deltap[deltap == 0.0] = np.inf
+        thermal = 1.0 / (np.exp(np.abs(deltap) / w_th) - 1.0)
+        thermal = [_data.Dense(thermal[:, :, k]) for k in range(2 * kmax + 1)]
+
+        gamma_kk = _data.add(gamma[0], gamma[0].transpose())
+        A = _data.add(gamma[0], _data.multiply(thermal[kmax], gamma_kk))
+
+        for k in range(1, kmax + 1):
+            g_kk = _data.add(gamma[k], gamma[-k].transpose())
+            thermal_kk = _data.multiply(thermal[kmax + k], g_kk)
+            A = _data.add(A, _data.add(gamma[k], thermal_kk))
+            thermal_kk = _data.multiply(thermal[kmax - k], g_kk.transpose())
+            A = _data.add(A, _data.add(gamma[-k], thermal_kk))
+    else:
+        # w_th is 0, thermal = 0s
+        A = gamma[0]
+        for k in range(1, kmax + 1):
+            A = _data.add(gamma[k], A)
+            A = _data.add(gamma[-k], A)
+
+    return A
+
+
+def _floquet_master_equation_tensor(A):
+    """
+    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
+    ----------
+    A : :class:`qutip.data.Data`
+        Floquet-Markov master equation rate matrix.
+
+    Returns
+    -------
+    output : array
+        The Floquet-Markov master equation tensor `R`.
+    """
+    N = A.shape[0]
+
+    # R[i+N*i, j+N*j] = A[j, i]
+    cols = np.arange(N, dtype=np.int32)
+    rows = np.linspace(1 - 1 / (N + 1), N, N**2 + 1, dtype=np.int32)
+    data = np.ones(N, dtype=complex)
+    expand = _data.csr.CSR((data, cols, rows), shape=(N**2, N))
+
+    R = expand @ A.transpose() @ expand.transpose()
+
+    # S[i+N*j, j+N*i] = -1/2 * sum_k(A[i, k] + A[j, k])
+    ket_1 = _data.Dense(np.ones(N, dtype=complex))
+    Asum = A @ ket_1
+    to_super = _data.add(_data.kron(Asum, ket_1), _data.kron(ket_1, Asum))
+    S = _data.diag(to_super.to_array().flatten() * -0.5, 0)
+
+    return _data.add(R, S)
+
+
+def floquet_tensor(H, c_ops, spectra_cb, T=0, w_th=0.0, kmax=5, nT=100):
+    """
+    Construct a tensor that represents the master equation in the floquet
+    basis.
+
+    Simplest RWA approximation [Grifoni et al, Phys.Rep. 304 229 (1998)]
+
+    Parameters
+    ----------
+    H : :class:`QobjEvo`
+        Periodic Hamiltonian
+
+    T : float
+        The period of the time-dependence of the hamiltonian.
+
+    c_ops : list of :class:`qutip.qobj`
+        list of collapse operators.
+
+    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`.
+
+    w_th : float
+        The temperature in units of frequency.
+
+    kmax : int
+        The truncation of the number of sidebands (default 5).
+
+    Returns
+    -------
+    output : array
+        The Floquet-Markov master equation tensor `R`.
+    """
+    if isinstance(H, FloquetBasis):
+        floquet_basis = H
+        T = H.T
+    else:
+        floquet_basis = FloquetBasis(H, T)
+    energy = floquet_basis.e_quasi
+    delta = _floquet_delta_tensor(energy, kmax, T)
+    x = _floquet_X_matrices(floquet_basis, c_ops, kmax, nT)
+    gamma = _floquet_gamma_matrices(x, delta, spectra_cb)
+    a = _floquet_A_matrix(delta, gamma, w_th)
+    r = _floquet_master_equation_tensor(a)
+    dims = floquet_basis.U(0).dims
+    return Qobj(
+        r, dims=[dims, dims], type="super", superrep="super", copy=False
+    )
+
+
+def fsesolve(H, psi0, tlist, e_ops=None, T=0.0, args=None, options=None):
+    """
+    Solve the Schrodinger equation using the Floquet formalism.
+
+    Parameters
+    ----------
+    H : :class:`Qobj`, :class:`QobjEvo`, :class:`QobjEvo` compatible format.
+        Periodic system Hamiltonian as :class:`QobjEvo`. List of
+        [:class:`Qobj`, :class:`Coefficient`] or callable that
+        can be made into :class:`QobjEvo` are also accepted.
+
+    psi0 : :class:`qutip.qobj`
+        Initial state vector (ket). If an operator is provided,
+
+    tlist : *list* / *array*
+        List of times for :math:`t`.
+
+    e_ops : list of :class:`qutip.qobj` / callback function, optional
+        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, default=tlist[-1]
+        The period of the time-dependence of the hamiltonian.
+
+    args : dictionary, optional
+        Dictionary with variables required to evaluate H.
+
+    options : dict, optional
+        Options for the results.
+
+        - store_final_state : bool
+          Whether or not to store the final state of the evolution in the
+          result class.
+        - store_states : bool, None
+          Whether or not to store the state vectors or density matrices.
+          On `None` the states will be saved if no expectation operators are
+          given.
+        - normalize_output : bool
+          Normalize output state to hide ODE numerical errors.
+
+    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 isinstance(H, FloquetBasis):
+        floquet_basis = H
+    else:
+        T = T or tlist[-1]
+        # `fsesolve` is a fallback from `fmmesolve`, for the later, options
+        # are for the open system evolution.
+        floquet_basis = FloquetBasis(H, T, args, precompute=tlist)
+
+    f_coeff = floquet_basis.to_floquet_basis(psi0)
+    result_options = {
+        "store_final_state": False,
+        "store_states": None,
+        "normalize_output": True,
+    }
+    result_options.update(options or {})
+    result = Result(e_ops, result_options, solver="fsesolve")
+    for t in tlist:
+        state_t = floquet_basis.from_floquet_basis(f_coeff, t)
+        result.add(t, state_t)
+
+    return result
+
+
+def fmmesolve(
+    H,
+    rho0,
+    tlist,
+    c_ops=None,
+    e_ops=None,
+    spectra_cb=None,
+    T=0,
+    w_th=0.0,
+    args=None,
+    options=None,
+):
+    """
+    Solve the dynamics for the system using the Floquet-Markov master equation.
+
+    Parameters
+    ----------
+    H : :class:`Qobj`, :class:`QobjEvo`, :class:`QobjEvo` compatible format.
+        Periodic system Hamiltonian as :class:`QobjEvo`. List of
+        [:class:`Qobj`, :class:`Coefficient`] or callable that
+        can be made into :class:`QobjEvo` are also accepted.
+
+    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. Time dependent collapse operators are not
+        supported.
+
+    e_ops : list of :class:`qutip.Qobj` / callback function
+        List of operators for which to evaluate expectation values.
+        The states are reverted to the lab basis before applying the
+
+    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.
+
+    w_th : float
+        The temperature of the environment in units of frequency.
+        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
+
+    args : *dictionary*
+        Dictionary of parameters for time-dependent Hamiltonian
+
+    options : None / dict
+        Dictionary of options for the solver.
+
+        - store_final_state : bool
+          Whether or not to store the final state of the evolution in the
+          result class.
+        - store_states : bool, None
+          Whether or not to store the state vectors or density matrices.
+          On `None` the states will be saved if no expectation operators are
+          given.
+        - store_floquet_states : bool
+          Whether or not to store the density matrices in the floquet basis in
+          ``result.floquet_states``.
+        - normalize_output : bool
+          Normalize output state to hide ODE numerical errors.
+        - progress_bar : str {'text', 'enhanced', 'tqdm', ''}
+          How to present the solver progress.
+          'tqdm' uses the python module of the same name and raise an error
+          if not installed. Empty string or False will disable the bar.
+        - progress_kwargs : dict
+          kwargs to pass to the progress_bar. Qutip's bars use `chunk_size`.
+        - method : str ["adams", "bdf", "lsoda", "dop853", "vern9", etc.]
+          Which differential equation integration method to use.
+        - atol, rtol : float
+          Absolute and relative tolerance of the ODE integrator.
+        - nsteps :
+          Maximum number of (internally defined) steps allowed in one ``tlist``
+          step.
+        - max_step : float, 0
+          Maximum lenght of one internal step. When using pulses, it should be
+          less than half the width of the thinnest pulse.
+
+        Other options could be supported depending on the integration method,
+        see `Integrator <./classes.html#classes-ode>`_.
+
+    Returns
+    -------
+    result: :class:`qutip.Result`
+
+        An instance of the class :class:`qutip.Result`, which contains
+        the expectation values for the times specified by `tlist`, and/or the
+        state density matrices corresponding to the times.
+    """
+    if c_ops is None:
+        return fsesolve(
+            H,
+            rho0,
+            tlist,
+            e_ops=e_ops,
+            T=T,
+            w_th=w_th,
+            args=args,
+            options=options,
+        )
+
+    if isinstance(H, FloquetBasis):
+        floquet_basis = H
+    else:
+        T = T or tlist[-1]
+        t_precompute = np.concatenate([tlist, np.linspace(0, T, 101)])
+        # `fsesolve` is a fallback from `fmmesolve`, for the later, options
+        # are for the open system evolution.
+        floquet_basis = FloquetBasis(H, T, args, precompute=t_precompute)
+
+    if not w_th and args:
+        w_th = args.get("w_th", 0.0)
+
+    if isinstance(c_ops, Qobj):
+        c_ops = [c_ops]
+
+    if spectra_cb is None:
+        spectra_cb = [lambda w: (w > 0)]
+    elif callable(spectra_cb):
+        spectra_cb = [spectra_cb]
+    if len(spectra_cb) == 1:
+        spectra_cb = spectra_cb * len(c_ops)
+
+    a_ops = list(zip(c_ops, spectra_cb))
+
+    solver = FMESolver(floquet_basis, a_ops, w_th=w_th, options=options)
+    return solver.run(rho0, tlist, e_ops=e_ops)
+
+
+class FloquetResult(Result):
+    def _post_init(self, floquet_basis):
+        self.floquet_basis = floquet_basis
+        if self.options["store_floquet_states"]:
+            self.floquet_states = []
+        else:
+            self.floquet_states = None
+        super()._post_init()
+
+    def add(self, t, state):
+        if self.options["store_floquet_states"]:
+            self.floquet_states.append(state)
+
+        state = self.floquet_basis.from_floquet_basis(state, t)
+        super().add(t, state)
+
+
+class FMESolver(MESolver):
+    """
+    Solver for the Floquet-Markov master equation.
+
+    .. note ::
+        Operators (``c_ops`` and ``e_ops``) are in the laboratory basis.
+
+    Parameters
+    ----------
+    floquet_basis : :class:`qutip.FloquetBasis`
+        The system Hamiltonian wrapped in a FloquetBasis object. Choosing a
+        different integrator for the ``floquet_basis`` than for the evolution
+        of the floquet state can improve the performance.
+
+    a_ops : list of tuple(:class:`qutip.Qobj`, callable)
+        List of collapse operators and the corresponding function for the noise
+        power spectrum. The collapse operator must be a :class:`Qobj` and
+        cannot be time dependent. The spectrum function must take and return
+        an numpy array.
+
+    w_th : float
+        The temperature of the environment in units of Hamiltonian frequency.
+
+    kmax : int [5]
+        The truncation of the number of sidebands..
+
+    nT : int [20*kmax]
+        The number of integration steps (for calculating X) within one period.
+
+    options : dict, optional
+        Options for the solver, see :obj:`FMESolver.options` and
+        `Integrator <./classes.html#classes-ode>`_ for a list of all options.
+    """
+
+    name = "fmmesolve"
+    _avail_integrators = {}
+    resultclass = FloquetResult
+    solver_options = {
+        "progress_bar": "text",
+        "progress_kwargs": {"chunk_size": 10},
+        "store_final_state": False,
+        "store_states": None,
+        "normalize_output": True,
+        "method": "adams",
+        "store_floquet_states": False,
+    }
+
+    def __init__(
+        self, floquet_basis, a_ops, w_th=0.0, *, kmax=5, nT=None, options=None
+    ):
+        self._options = {}
+        self.options = {} if options is None else options
+        if isinstance(floquet_basis, FloquetBasis):
+            self.floquet_basis = floquet_basis
+        else:
+            raise TypeError("The ``floquet_basis`` must be a FloquetBasis")
+
+        nT = nT or max(100, 20 * kmax)
+        self._num_collapse = len(a_ops)
+        c_ops, spectra_cb = zip(*a_ops)
+        if not all(
+            isinstance(c_op, Qobj) and callable(spectrum)
+            for c_op, spectrum in a_ops
+        ):
+            raise TypeError("a_ops must be tuple of (Qobj, callable)")
+        self.rhs = QobjEvo(
+            floquet_tensor(
+                self.floquet_basis,
+                c_ops,
+                spectra_cb,
+                w_th=w_th,
+                kmax=kmax,
+                nT=nT,
+            )
+        )
+
+        self._integrator = self._get_integrator()
+        self._state_metadata = {}
+        self.stats = self._initialize_stats()
+
+    def _initialize_stats(self):
+        stats = Solver._initialize_stats(self)
+        stats.update(
+            {
+                "solver": "Floquet-Markov master equation",
+                "num_collapse": self._num_collapse,
+            }
+        )
+        return stats
+
+    def _argument(self, args):
+        if args:
+            raise ValueError("FMESolver cannot update arguments")
+
+    def start(self, state0, t0, *, floquet=False):
+        """
+        Set the initial state and time for a step evolution.
+        ``options`` for the evolutions are read at this step.
+
+        Parameters
+        ----------
+        state0 : :class:`Qobj`
+            Initial state of the evolution.
+
+        t0 : double
+            Initial time of the evolution.
+
+        floquet : bool, optional {False}
+            Whether the initial state is in the floquet basis or laboratory
+            basis.
+        """
+        if not floquet:
+            state0 = self.floquet_basis.to_floquet_basis(state0, t0)
+        super().start(state0, t0)
+
+    def step(self, t, *, args=None, copy=True, floquet=False):
+        """
+        Evolve the state to ``t`` and return the state as a :class:`Qobj`.
+
+        Parameters
+        ----------
+        t : double
+            Time to evolve to, must be higher than the last call.
+
+        copy : bool, optional {True}
+            Whether to return a copy of the data or the data in the ODE solver.
+
+        floquet : bool, optional {False}
+            Whether to return the state in the floquet basis or laboratory
+            basis.
+
+        args : dict, optional {None}
+            Not supported
+
+        .. note::
+            The state must be initialized first by calling ``start`` or
+            ``run``. If ``run`` is called, ``step`` will continue from the last
+            time and state obtained.
+        """
+        if args:
+            raise ValueError("FMESolver cannot update arguments")
+        state = super().step(t)
+        if not floquet:
+            state = self.floquet_basis.from_floquet_basis(state, t)
+        elif copy:
+            state = state.copy()
+        return state
+
+    def run(self, state0, tlist, *, floquet=False, args=None, e_ops=None):
+        """
+        Calculate the evolution of the quantum system.
+
+        For a ``state0`` at time ``tlist[0]`` do the evolution as directed by
+        ``rhs`` and for each time in ``tlist`` store the state and/or
+        expectation values in a :class:`Result`. The evolution method and
+        stored results are determined by ``options``.
+
+        Parameters
+        ----------
+        state0 : :class:`Qobj`
+            Initial state of the evolution.
+
+        tlist : list of double
+            Time for which to save the results (state and/or expect) of the
+            evolution. The first element of the list is the initial time of the
+            evolution. Each times of the list must be increasing, but does not
+            need to be uniformy distributed.
+
+        floquet : bool, optional {False}
+            Whether the initial state in the floquet basis or laboratory basis.
+
+        args : dict, optional {None}
+            Not supported
+
+        e_ops : list {None}
+            List of Qobj, QobjEvo or callable to compute the expectation
+            values. Function[s] must have the signature
+            f(t : float, state : Qobj) -> expect.
+
+        Returns
+        -------
+        results : :class:`qutip.solver.FloquetResult`
+            Results of the evolution. States and/or expect will be saved. You
+            can control the saved data in the options.
+        """
+        if args:
+            raise ValueError("FMESolver cannot update arguments")
+        if not floquet:
+            state0 = self.floquet_basis.to_floquet_basis(state0, tlist[0])
+        _time_start = time()
+        _data0 = self._prepare_state(state0)
+        self._integrator.set_state(tlist[0], _data0)
+        stats = self._initialize_stats()
+        results = self.resultclass(
+            e_ops,
+            self.options,
+            solver=self.name,
+            stats=stats,
+            floquet_basis=self.floquet_basis,
+        )
+        results.add(tlist[0], self._restore_state(_data0, copy=False))
+        stats["preparation time"] += time() - _time_start
+
+        progress_bar = progess_bars[self.options["progress_bar"]]()
+        progress_bar.start(len(tlist) - 1, **self.options["progress_kwargs"])
+        for t, state in self._integrator.run(tlist):
+            progress_bar.update()
+            results.add(t, self._restore_state(state, copy=False))
+        progress_bar.finished()
+
+        stats["run time"] = progress_bar.total_time()
+        # TODO: It would be nice if integrator could give evolution statistics
+        # stats.update(_integrator.stats)
+        return results
diff --git a/qutip/solver/floquet_bwcomp.py b/qutip/solver/floquet_bwcomp.py
new file mode 100644
index 0000000000..05a5fce871
--- /dev/null
+++ b/qutip/solver/floquet_bwcomp.py
@@ -0,0 +1,230 @@
+""" Floquet solver compatibility functions that behave like the corresponding
+functions from QuTiP 4.7.
+
+These functions are indented to be used when porting code from QuTiP 4.7 to
+QuTiP 5. They are deprecated and will be removed in QuTiP 5.1.
+"""
+
+__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",
+    "floquet_master_equation_rates",
+]
+
+from .floquet import *
+import numpy as np
+import warnings
+
+
+def floquet_modes(H, T, args=None, sort=False, U=None, options=None):
+    """
+    Calculate the initial Floquet modes Phi_alpha(0) for a driven system with
+    period T.
+
+    Deprecated from qutip v5. Use :class:`FloquetBasis` instead:
+
+    fbasis = FloquetBasis(H, T, args=args, options=options, sort=sort)
+    f_mode_0 = fbasis.mode(0)
+    f_energies = fbasis.e_quasi
+    """
+    warnings.warn(FutureWarning("`floquet_modes` is deprecated."))
+    fbasis = FloquetBasis(H, T, args=args, options=options, sort=sort)
+    f_mode_0 = fbasis.mode(0)
+    f_energies = fbasis.e_quasi
+    return f_mode_0, f_energies
+
+
+def floquet_modes_t(f_modes_0, f_energies, t, H, T, args=None, options=None):
+    """
+    Calculate the Floquet modes at times tlist Phi_alpha(tlist) propagting the
+    initial Floquet modes Phi_alpha(0).
+
+    Deprecated from qutip v5. Use :class:`FloquetBasis` instead:
+
+    fbasis = FloquetBasis(H, T, args=args, options=options)
+    f_mode_t = fbasis.mode(t)
+    """
+    warnings.warn(FutureWarning("`floquet_modes_t` is deprecated."))
+    fbasis = FloquetBasis(H, T, args=args, options=options)
+    return fbasis.mode(t)
+
+
+def floquet_modes_table(
+    f_modes_0, f_energies, tlist, H, T, args=None, options=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.
+
+    Deprecated from qutip v5. Use :class:`FloquetBasis` instead:
+
+    fbasis = FloquetBasis(H, T, args=args, options=options, precompute=tlist)
+    """
+    warnings.warn(FutureWarning("`floquet_modes_table` is deprecated."))
+    return FloquetBasis(H, T, args=args, options=options, precompute=tlist)
+
+
+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.
+
+    Deprecated from qutip v5. Use :class:`FloquetBasis` instead:
+
+    f_modes_table_t = fbasis = FloquetBasis(...)
+    f_mode_t = f_modes_table_t.mode(t)
+    """
+    warnings.warn(FutureWarning("`floquet_modes_t_lookup` is deprecated."))
+    return f_modes_table_t.mode(t)
+
+
+def floquet_states(f_modes_t, f_energies, t):
+    """
+    Evaluate the floquet states at time t given the Floquet modes at that time.
+
+    Deprecated from qutip v5. Use :class:`FloquetBasis` instead:
+
+    fbasis = FloquetBasis(H, T, args=args, options=options)
+    f_state_t = fbasis.state(t)
+    """
+    warnings.warn(FutureWarning("`floquet_states` is deprecated."))
+    return [
+        (f_modes_t[i] * np.exp(-1j * f_energies[i] * t))
+        for i in np.arange(len(f_energies))
+    ]
+
+
+def floquet_states_t(f_modes_0, f_energies, t, H, T, args=None, options=None):
+    """
+    Evaluate the floquet states at time t given the initial Floquet modes.
+
+    Deprecated from qutip v5. Use :class:`FloquetBasis` instead:
+
+    fbasis = FloquetBasis(H, T, args=args, options=options)
+    f_state_t = fbasis.state(t)
+    """
+    warnings.warn(FutureWarning("`floquet_states` is deprecated."))
+    fbasis = FloquetBasis(H, T, args=args, options=options)
+    return fbasis.state(t)
+
+
+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`.
+
+    Deprecated from qutip v5. Use :class:`FloquetBasis` instead:
+
+    fbasis = FloquetBasis(H, T, args=args, options=options)
+    psi_t = fbasis.from_floquet_basis(f_coeff, t)
+    """
+    warnings.warn(FutureWarning("`floquet_wavefunction` is deprecated."))
+    return sum(
+        [
+            f_modes_t[i] * np.exp(-1j * f_energies[i] * t) * f_coeff[i]
+            for i in np.arange(1, len(f_energies))
+        ],
+        start=f_modes_t[0] * np.exp(-1j * f_energies[0] * t) * f_coeff[0],
+    )
+
+
+def floquet_wavefunction_t(
+    f_modes_0, f_energies, f_coeff, t, H, T, args=None, options=None
+):
+    """
+    Evaluate the wavefunction for a time t using the Floquet state
+    decompositon, given the initial Floquet modes.
+
+    Deprecated from qutip v5. Use :class:`FloquetBasis` instead:
+
+    fbasis = FloquetBasis(H, T, args=args, options=options)
+    psi_t = fbasis.from_floquet_basis(f_coeff, t)
+    """
+    warnings.warn(FutureWarning("`floquet_states` is deprecated."))
+    fbasis = FloquetBasis(H, T, args=args, options=options)
+    return fbasis.from_floquet_basis(f_coeff, t)
+
+
+def floquet_state_decomposition(f_states, f_energies, psi):
+    r"""
+    Decompose the wavefunction `psi` (typically an initial state) in terms of
+    the Floquet states, :math:`\psi = \sum_\alpha c_\alpha \psi_\alpha(0)`.
+
+    Deprecated from qutip v5. Use :class:`FloquetBasis` instead:
+
+    fbasis = FloquetBasis(H, T, args=args, options=options)
+    f_coeff = fbasis.to_floquet_basis(psi)
+    """
+    warnings.warn(FutureWarning("`floquet_states` is deprecated."))
+    return [state.dag() * psi for state in f_states]
+
+
+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.
+
+    .. note ::
+
+        Deprecated. For the Floquet-Markov master equation's tensor, use
+        :func:`floquet_tensor`. For the rates matrices, use
+        :func:`floquet_delta_tensor`, :func:`floquet_X_matrices`,
+        :func:`floquet_gamma_matrices` and/or s:func:`floquet_A_matrix`.
+
+    Parameters
+    ----------
+    f_modes_0 : Any
+        No longer used.
+    f_energies : Any
+        No longer used.
+    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.
+    kmax : int, default=5
+        The truncation of the number of sidebands.
+    f_modes_table_t : Any
+        No longer used.
+
+    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.
+    """
+    warnings.warn(FutureWarning("`floquet_master_equation_rates` is deprecated."))
+    floquet_basis = FloquetBasis(H, T, args=args)
+    energy = floquet_basis.e_quasi
+    delta = floquet_delta_tensor(energy, kmax, T)
+    x = floquet_X_matrices(floquet_basis, [c_op], kmax, nT)
+    gamma = floquet_gamma_matrices(x, delta, [J_cb])
+    a = floquet_A_matrix(delta, gamma, w_th)
+    return delta, x[0], gamma, a
diff --git a/qutip/solver/propagator.py b/qutip/solver/propagator.py
index fbc83698ad..179c0dcef0 100644
--- a/qutip/solver/propagator.py
+++ b/qutip/solver/propagator.py
@@ -169,6 +169,7 @@ def __init__(self, system, *, c_ops=(), args=None, options=None,
         self.times = [0]
         self.invs = [None]
         self.props = [qeye(self.solver.sys_dims)]
+        self.solver.start(self.props[0], self.times[0])
         self.cte = self.solver.rhs.isconstant
         H_0 = self.solver.rhs(0)
         self.unitary = not H_0.issuper and H_0.isherm
@@ -209,8 +210,10 @@ def __call__(self, t, t_start=0, **args):
         # We could improve it when the system is constant using U(2t) = U(t)**2
         if not self.cte and args and args != self.args:
             self.args = args
+            self.solver._argument(args)
             self.times = [0]
             self.props = [qeye(self.props[0].dims[0])]
+            self.solver.start(self.props[0], self.times[0])
 
         if t_start:
             if t == t_start:
@@ -245,13 +248,16 @@ def _compute(self, t, idx):
         Compute the propagator at ``t``, ``idx`` point to a pair of
         (time, propagator) close to the desired time.
         """
-        if idx > 0:
+        t_last = self.solver._integrator.get_state(copy=False)[0]
+        if self.times[idx-1] <= t_last <= t:
+            U = self.solver.step(t)
+        elif idx > 0:
             self.solver.start(self.props[idx-1], self.times[idx-1])
-            U = self.solver.step(t, args=self.args)
+            U = self.solver.step(t)
         else:
             # Evolving backward in time is not supported by all integrator.
             self.solver.start(qeye(self.props[0].dims[0]), t)
-            Uinv = self.solver.step(self.times[idx], args=self.args)
+            Uinv = self.solver.step(self.times[idx])
             U = self._inv(Uinv)
         return U
 
diff --git a/qutip/tests/solve/test_floquet.py b/qutip/tests/solve/test_floquet.py
deleted file mode 100644
index f4b8c09d84..0000000000
--- a/qutip/tests/solve/test_floquet.py
+++ /dev/null
@@ -1,413 +0,0 @@
-import numpy as np
-from qutip import fsesolve, sigmax, sigmaz, rand_ket, num, mesolve, sigmay
-from qutip import sigmap, sigmam, floquet_master_equation_rates, expect, Qobj
-from qutip import floquet_modes, floquet_modes_table, fmmesolve
-from qutip import floquet_modes_t_lookup
-import pytest
-
-
-class TestFloquet:
-    """
-    A test class for the QuTiP functions for Floquet formalism.
-    """
-
-    def testFloquetUnitary(self):
-        """
-        Floquet: test unitary evolution of time-dependent two-level system
-        """
-
-        delta = 1.0 * 2 * np.pi
-        eps0 = 1.0 * 2 * np.pi
-        A = 0.5 * 2 * np.pi
-        omega = np.sqrt(delta ** 2 + eps0 ** 2)
-        T = (2 * np.pi) / omega
-        tlist = np.linspace(0.0, 2 * T, 101)
-        psi0 = rand_ket(2)
-        H0 = - eps0 / 2.0 * sigmaz() - delta / 2.0 * sigmax()
-        H1 = A / 2.0 * sigmax()
-        args = {'w': omega}
-        H = [H0, [H1, lambda t, args: np.sin(args['w'] * t)]]
-        e_ops = [num(2)]
-
-        # Solve schrodinger equation with floquet solver
-        sol = fsesolve(H, psi0, tlist, e_ops, T, args)
-
-        # Compare with results from standard schrodinger equation
-        sol_ref = mesolve(H, psi0, tlist, [], e_ops, args)
-
-        np.testing.assert_allclose(sol.expect[0], sol_ref.expect[0], atol=1e-4)
-
-    def testFloquetMasterEquation1(self):
-        """
-        Test Floquet-Markov Master Equation for a driven two-level system
-        without dissipation.
-        """
-
-        delta = 1.0 * 2 * np.pi
-        eps0 = 1.0 * 2 * np.pi
-        A = 0.5 * 2 * np.pi
-        omega = np.sqrt(delta ** 2 + eps0 ** 2)
-        T = (2 * np.pi) / omega
-        tlist = np.linspace(0.0, 2 * T, 101)
-        psi0 = rand_ket(2)
-        H0 = - eps0 / 2.0 * sigmaz() - delta / 2.0 * sigmax()
-        H1 = A / 2.0 * sigmax()
-        args = {'w': omega}
-        H = [H0, [H1, lambda t, args: np.sin(args['w'] * t)]]
-        e_ops = [num(2)]
-        gamma1 = 0
-
-        # Collapse operator for Floquet-Markov Master Equation
-        c_op_fmmesolve = sigmax()
-
-        # Collapse operators for Lindblad Master Equation
-        def noise_spectrum(omega):
-            if omega > 0:
-                return 0.5 * gamma1 * omega/(2*np.pi)
-            else:
-                return 0
-
-        ep, vp = H0.eigenstates()
-        op0 = vp[0]*vp[0].dag()
-        op1 = vp[1]*vp[1].dag()
-
-        c_op_mesolve = []
-        for i in range(2):
-            for j in range(2):
-                if i != j:
-                    # caclculate the rate
-                    gamma = 2 * np.pi * c_op_fmmesolve.matrix_element(
-                        vp[j], vp[i]) * c_op_fmmesolve.matrix_element(
-                        vp[i], vp[j]) * noise_spectrum(ep[j] - ep[i])
-
-                    # add c_op for mesolve
-                    c_op_mesolve.append(
-                        np.sqrt(gamma) * (vp[i] * vp[j].dag())
-                    )
-
-        # Find the floquet modes
-        f_modes_0, f_energies = floquet_modes(H, T, args)
-
-        # Precalculate mode table
-        f_modes_table_t = floquet_modes_table(f_modes_0, f_energies,
-                                              np.linspace(0, T, 500 + 1),
-                                              H, T, args)
-
-        # Solve the floquet-markov master equation
-        output1 = fmmesolve(H, psi0, tlist, [c_op_fmmesolve], [],
-                            [noise_spectrum], T, args, floquet_basis=True)
-
-        # Calculate expectation values in the computational basis
-        p_ex = np.zeros(np.shape(tlist), dtype=complex)
-        for idx, t in enumerate(tlist):
-            f_modes_t = floquet_modes_t_lookup(f_modes_table_t, t, T)
-            f_states_t = [np.exp(-1j*t*f_energies[0])*f_modes_t[0],
-                          np.exp(-1j*t*f_energies[1])*f_modes_t[1]]
-            p_ex[idx] = expect(num(2), output1.states[idx].transform(
-                f_states_t, True))
-
-        # Compare with mesolve
-        output2 = mesolve(H, psi0, tlist, c_op_mesolve, [], args)
-        p_ex_ref = expect(num(2), output2.states)
-
-        np.testing.assert_allclose(np.real(p_ex), np.real(p_ex_ref), atol=1e-4)
-
-    @pytest.mark.parametrize("kmax", [5, 100, 200])
-    def testFloquetMasterEquation2(self, kmax):
-        """
-        Test Floquet-Markov Master Equation for a two-level system
-        subject to dissipation.
-        """
-
-        delta = 1.0 * 2 * np.pi
-        eps0 = 1.0 * 2 * np.pi
-        A = 0.5 * 2 * np.pi
-        omega = np.sqrt(delta ** 2 + eps0 ** 2)
-        T = (2 * np.pi) / omega
-        tlist = np.linspace(0.0, 2 * T, 101)
-        psi0 = rand_ket(2)
-        H0 = - eps0 / 2.0 * sigmaz() - delta / 2.0 * sigmax()
-        H1 = A / 2.0 * sigmax()
-        args = {'w': omega}
-        H = [H0, [H1, lambda t, args: np.sin(args['w'] * t)]]
-        e_ops = [num(2)]
-        gamma1 = 1
-
-        A = 0. * 2 * np.pi
-        psi0 = rand_ket(2)
-        H1 = A / 2.0 * sigmax()
-        args = {'w': omega}
-        H = [H0, [H1, lambda t, args: np.sin(args['w'] * t)]]
-
-        # Collapse operator for Floquet-Markov Master Equation
-        c_op_fmmesolve = sigmax()
-
-        # Collapse operator for Lindblad Master Equation
-        def noise_spectrum(omega):
-            if omega > 0:
-                return 0.5 * gamma1 * omega/(2*np.pi)
-            else:
-                return 0
-
-        ep, vp = H0.eigenstates()
-        op0 = vp[0]*vp[0].dag()
-        op1 = vp[1]*vp[1].dag()
-
-        c_op_mesolve = []
-        for i in range(2):
-            for j in range(2):
-                if i != j:
-                    # caclculate the rate
-                    gamma = 2 * np.pi * c_op_fmmesolve.matrix_element(
-                        vp[j], vp[i]) * c_op_fmmesolve.matrix_element(
-                        vp[i], vp[j]) * noise_spectrum(ep[j] - ep[i])
-
-                    # add c_op for mesolve
-                    c_op_mesolve.append(
-                        np.sqrt(gamma) * (vp[i] * vp[j].dag())
-                    )
-
-        # Find the floquet modes
-        f_modes_0, f_energies = floquet_modes(H, T, args)
-
-        # Precalculate mode table
-        f_modes_table_t = floquet_modes_table(f_modes_0, f_energies,
-                                              np.linspace(0, T, 500 + 1),
-                                              H, T, args)
-
-        # Solve the floquet-markov master equation
-        output1 = fmmesolve(
-                            H, psi0, tlist, [c_op_fmmesolve], [],
-                            [noise_spectrum], T, args, floquet_basis=True,
-                            kmax=kmax)
-        # Calculate expectation values in the computational basis
-        p_ex = np.zeros(np.shape(tlist), dtype=complex)
-        for idx, t in enumerate(tlist):
-            f_modes_t = floquet_modes_t_lookup(f_modes_table_t, t, T)
-            f_states_t = [np.exp(-1j*t*f_energies[0])*f_modes_t[0],
-                          np.exp(-1j*t*f_energies[1])*f_modes_t[1]]
-            p_ex[idx] = expect(num(2), output1.states[idx].transform(
-                                f_states_t,
-                                True))
-
-        # Compare with mesolve
-        output2 = mesolve(H, psi0, tlist, c_op_mesolve, [], args)
-        p_ex_ref = expect(num(2), output2.states)
-
-        np.testing.assert_allclose(np.real(p_ex), np.real(p_ex_ref), atol=1e-4)
-
-    @pytest.mark.parametrize("kmax", [5, 100, 300])
-    def testFloquetMasterEquation3(self, kmax):
-        """
-        Test Floquet-Markov Master Equation for a two-level system
-        subject to dissipation with internal transform of fmmesolve
-        """
-
-        delta = 1.0 * 2 * np.pi
-        eps0 = 1.0 * 2 * np.pi
-        A = 0.5 * 2 * np.pi
-        omega = np.sqrt(delta ** 2 + eps0 ** 2)
-        T = (2 * np.pi) / omega
-        tlist = np.linspace(0.0, 2 * T, 101)
-        psi0 = rand_ket(2)
-        H0 = - eps0 / 2.0 * sigmaz() - delta / 2.0 * sigmax()
-        H1 = A / 2.0 * sigmax()
-        args = {'w': omega}
-        H = [H0, [H1, lambda t, args: np.sin(args['w'] * t)]]
-        e_ops = [num(2)]
-        gamma1 = 1
-
-        A = 0. * 2 * np.pi
-        psi0 = rand_ket(2)
-        H1 = A / 2.0 * sigmax()
-        args = {'w': omega}
-        H = [H0, [H1, lambda t, args: np.sin(args['w'] * t)]]
-
-        # Collapse operator for Floquet-Markov Master Equation
-        c_op_fmmesolve = sigmax()
-
-        # Collapse operator for Lindblad Master Equation
-        def noise_spectrum(omega):
-            if omega > 0:
-                return 0.5 * gamma1 * omega/(2*np.pi)
-            else:
-                return 0
-
-        ep, vp = H0.eigenstates()
-        op0 = vp[0]*vp[0].dag()
-        op1 = vp[1]*vp[1].dag()
-
-        c_op_mesolve = []
-        for i in range(2):
-            for j in range(2):
-                if i != j:
-                    # caclculate the rate
-                    gamma = 2*np.pi*c_op_fmmesolve.matrix_element(
-                        vp[j], vp[i])*c_op_fmmesolve.matrix_element(
-                        vp[i], vp[j])*noise_spectrum(ep[j]-ep[i])
-
-                    # add c_op for mesolve
-                    c_op_mesolve.append(
-                        np.sqrt(gamma) * (vp[i] * vp[j].dag())
-                    )
-
-        # Solve the floquet-markov master equation
-        output1 = fmmesolve(
-                            H, psi0, tlist, [c_op_fmmesolve], [num(2)],
-                            [noise_spectrum], T, args, floquet_basis=False,
-                            kmax=kmax)
-        p_ex = output1.expect[0]
-        # Compare with mesolve
-        output2 = mesolve(H, psi0, tlist, c_op_mesolve, [num(2)], args)
-        p_ex_ref = output2.expect[0]
-
-        np.testing.assert_allclose(np.real(p_ex), np.real(p_ex_ref),
-                                   atol=5 * 1e-4)
-
-    def testFloquetMasterEquation_multiple_coupling(self):
-        """
-        Test Floquet-Markov Master Equation for a two-level system
-        subject to dissipation with multiple coupling operators
-        """
-
-        delta = 1.0 * 2 * np.pi
-        eps0 = 1.0 * 2 * np.pi
-        A = 0.5 * 2 * np.pi
-        omega = np.sqrt(delta ** 2 + eps0 ** 2)
-        T = (2 * np.pi) / omega
-        tlist = np.linspace(0.0, 2 * T, 101)
-        psi0 = rand_ket(2)
-        H0 = - eps0 / 2.0 * sigmaz() - delta / 2.0 * sigmax()
-        H1 = A / 2.0 * sigmax()
-        args = {'w': omega}
-        H = [H0, [H1, lambda t, args: np.sin(args['w'] * t)]]
-        e_ops = [num(2)]
-        gamma1 = 1
-
-        A = 0. * 2 * np.pi
-        psi0 = rand_ket(2)
-        H1 = A / 2.0 * sigmax()
-        args = {'w': omega}
-        H = [H0, [H1, lambda t, args: np.sin(args['w'] * t)]]
-
-        # Collapse operator for Floquet-Markov Master Equation
-        c_ops_fmmesolve = [sigmax(), sigmay()]
-
-        # Collapse operator for Lindblad Master Equation
-        def noise_spectrum1(omega):
-            if omega > 0:
-                return 0.5 * gamma1 * omega/(2*np.pi)
-            else:
-                return 0
-        def noise_spectrum2(omega):
-            if omega > 0:
-                return 0.5 * gamma1 / (2 * np.pi)
-            else:
-                return 0
-
-        noise_spectra = [noise_spectrum1, noise_spectrum2]
-
-        ep, vp = H0.eigenstates()
-        op0 = vp[0]*vp[0].dag()
-        op1 = vp[1]*vp[1].dag()
-
-        c_op_mesolve = []
-
-        # Convert the c_ops for fmmesolve to c_ops for mesolve
-        for c_op_fmmesolve, noise_spectrum in zip(c_ops_fmmesolve,
-                                                   noise_spectra):
-            for i in range(2):
-                for j in range(2):
-                    if i != j:
-                        # caclculate the rate
-                        gamma = 2*np.pi*c_op_fmmesolve.matrix_element(
-                            vp[j], vp[i])*c_op_fmmesolve.matrix_element(
-                            vp[i], vp[j])*noise_spectrum(ep[j]-ep[i])
-                        # add c_op for mesolve
-                        c_op_mesolve.append(
-                            np.sqrt(gamma) * (vp[i] * vp[j].dag())
-                        )
-
-        # Solve the floquet-markov master equation
-        output1 = fmmesolve(
-                            H, psi0, tlist, c_ops_fmmesolve, [num(2)],
-                            noise_spectra, T, args, floquet_basis=False)
-        p_ex = output1.expect[0]
-        # Compare with mesolve
-        output2 = mesolve(H, psi0, tlist, c_op_mesolve, [num(2)], args)
-        p_ex_ref = output2.expect[0]
-
-        np.testing.assert_allclose(np.real(p_ex), np.real(p_ex_ref), atol=1e-4)
-
-    def testFloquetRates(self):
-        """
-        Compare rate transition and frequency transitions to analytical
-        results for a driven two-level system, for different drive amplitudes.
-        """
-
-        # Parameters
-        wq = 5.4 * 2 * np.pi
-        wd = 6.7 * 2 * np.pi
-        delta = wq - wd
-        T = 2 * np.pi / wd
-        tlist = np.linspace(0.0, 2 * T, 101)
-        array_A = np.linspace(0.001, 3 * 2 * np.pi, 10, endpoint=False)
-        H0 = Qobj(wq/2 * sigmaz())
-        arg = {'wd': wd}
-        c_ops = sigmax()
-        gamma1 = 1
-
-        delta_ana_deltas = []
-        array_ana_E0 = [-np.sqrt((delta/2)**2 + a**2) for a in array_A]
-        array_ana_E1 = [np.sqrt((delta/2)**2 + a**2) for a in array_A]
-        array_ana_delta = [2*np.sqrt((delta/2)**2 + a**2) for a in array_A]
-
-        def noise_spectrum(omega):
-            if omega > 0:
-                return 0.5 * gamma1 * omega/(2*np.pi)
-            else:
-                return 0
-
-        idx = 0
-        for a in array_A:
-            # Hamiltonian
-            H1_p = Qobj(a * sigmap())
-            H1_m = Qobj(a * sigmam())
-            H = [H0, [H1_p, lambda t, args: np.exp(-1j * arg['wd'] * t)],
-                     [H1_m, lambda t, args: np.exp(1j * arg['wd'] * t)]]
-
-            # Floquet modes
-            fmodes0, fenergies = floquet_modes(H, T, args={}, sort=True)
-            f_modes_table_t = floquet_modes_table(fmodes0, fenergies,
-                                                  tlist, H, T,
-                                                  args={})
-            # Get X delta
-            DeltaMatrix, X, frates, Amat = floquet_master_equation_rates(
-                                            fmodes0, fenergies,
-                                            c_ops, H, T, {},
-                                            noise_spectrum, 0, 5)
-            # Check energies
-            deltas = np.ndarray.flatten(DeltaMatrix)
-
-            # deltas and array_ana_delta have at least 1 value in common.
-            assert (min(abs(deltas-array_ana_delta[idx])) < 1e-4)
-
-            # Check matrix elements
-            Xs = np.ndarray.flatten(X)
-
-            normPlus = np.sqrt(a**2 + (array_ana_E1[idx] - delta/2)**2)
-            normMinus = np.sqrt(a**2 + (array_ana_E0[idx] - delta/2)**2)
-
-            Xpp_p1 = (a/normPlus**2)*(array_ana_E1[idx]-delta/2)
-            assert (min(abs(Xs-Xpp_p1)) < 1e-4)
-            Xpp_m1 = (a/normPlus**2)*(array_ana_E1[idx]-delta/2)
-            assert (min(abs(Xs-Xpp_m1)) < 1e-4)
-            Xmm_p1 = (a/normMinus**2)*(array_ana_E0[idx]-delta/2)
-            assert (min(abs(Xs-Xmm_p1)) < 1e-4)
-            Xmm_m1 = (a/normMinus**2)*(array_ana_E0[idx]-delta/2)
-            assert (min(abs(Xs-Xmm_m1)) < 1e-4)
-            Xpm_p1 = (a/(normMinus*normPlus))*(array_ana_E0[idx]-delta/2)
-            assert (min(abs(Xs-Xmm_p1)) < 1e-4)
-            Xpm_m1 = (a/(normMinus*normPlus))*(array_ana_E1[idx]-delta/2)
-            assert (min(abs(Xs-Xpm_m1)) < 1e-4)
-            idx += 1
diff --git a/qutip/tests/solver/test_floquet.py b/qutip/tests/solver/test_floquet.py
new file mode 100644
index 0000000000..be90e6fab2
--- /dev/null
+++ b/qutip/tests/solver/test_floquet.py
@@ -0,0 +1,333 @@
+import numpy as np
+from qutip import (
+    sigmax, sigmay, sigmaz,  sigmap, sigmam,
+    rand_ket, num, destroy,
+    mesolve, expect, sesolve,
+    Qobj, QobjEvo, coefficient
+ )
+
+from qutip.solver.floquet import (
+    FloquetBasis, floquet_tensor, fmmesolve, FMESolver,
+    _floquet_delta_tensor, _floquet_X_matrices, fsesolve
+)
+import pytest
+
+
+def _convert_c_ops(c_op_fmmesolve, noise_spectrum, vp, ep):
+    """
+    Convert e_ops for fmmesolve to mesolve
+    """
+    c_op_mesolve = []
+    N = len(vp)
+    for i in range(N):
+        for j in range(N):
+            if i != j:
+                # caclculate the rate
+                gamma = 2 * np.pi * c_op_fmmesolve.matrix_element(
+                    vp[j], vp[i]) * c_op_fmmesolve.matrix_element(
+                    vp[i], vp[j]) * noise_spectrum(ep[j] - ep[i])
+
+                # add c_op for mesolve
+                c_op_mesolve.append(
+                    np.sqrt(gamma) * (vp[i] * vp[j].dag())
+                )
+    return c_op_mesolve
+
+
+class TestFloquet:
+    """
+    A test class for the QuTiP functions for Floquet formalism.
+    """
+
+    def testFloquetBasis(self):
+        N = 10
+        a = destroy(N)
+        H = num(N) + (a+a.dag()) * coefficient(lambda t: np.cos(t))
+        T = 2 * np.pi
+        floquet_basis = FloquetBasis(H, T)
+        psi0 = rand_ket(N)
+        tlist = np.linspace(0, 10, 11)
+        floquet_psi0 = floquet_basis.to_floquet_basis(psi0)
+        states = sesolve(H, psi0, tlist).states
+        for t, state in zip(tlist, states):
+            from_floquet = floquet_basis.from_floquet_basis(floquet_psi0, t)
+            assert state.overlap(from_floquet) == pytest.approx(1., abs=5e-5)
+
+    def testFloquetUnitary(self):
+        N = 10
+        a = destroy(N)
+        H = num(N) + (a+a.dag()) * coefficient(lambda t: np.cos(t))
+        T = 2 * np.pi
+        psi0 = rand_ket(N)
+        tlist = np.linspace(0, 10, 11)
+        states_se = sesolve(H, psi0, tlist).states
+        states_fse = fsesolve(H, psi0, tlist, T=T).states
+        for state_se, state_fse in zip(states_se, states_fse):
+            assert state_se.overlap(state_fse) == pytest.approx(1., abs=5e-5)
+
+
+    def testFloquetMasterEquation1(self):
+        """
+        Test Floquet-Markov Master Equation for a driven two-level system
+        without dissipation.
+        """
+        delta = 1.0 * 2 * np.pi
+        eps0 = 1.0 * 2 * np.pi
+        A = 0.5 * 2 * np.pi
+        omega = np.sqrt(delta ** 2 + eps0 ** 2)
+        T = (2 * np.pi) / omega
+        tlist = np.linspace(0.0, 2 * T, 101)
+        psi0 = rand_ket(2)
+        H0 = - eps0 / 2.0 * sigmaz() - delta / 2.0 * sigmax()
+        H1 = A / 2.0 * sigmax()
+        args = {'w': omega}
+        H = [H0, [H1, lambda t, args: np.sin(args['w'] * t)]]
+        e_ops = [num(2)]
+        gamma1 = 0
+
+        # Collapse operator for Floquet-Markov Master Equation
+        c_op_fmmesolve = sigmax()
+
+        # Collapse operators for Lindblad Master Equation
+        def spectrum(omega):
+            return (omega > 0) * omega * 0.5 * gamma1 / (2 * np.pi)
+
+        ep, vp = H0.eigenstates()
+        op0 = vp[0]*vp[0].dag()
+        op1 = vp[1]*vp[1].dag()
+
+        c_op_mesolve = _convert_c_ops(c_op_fmmesolve, spectrum, vp, ep)
+
+        # Solve the floquet-markov master equation
+        p_ex = fmmesolve(H, psi0, tlist, [c_op_fmmesolve], [num(2)],
+                         [spectrum], T, args=args).expect[0]
+
+        # Compare with mesolve
+        p_ex_ref = mesolve(
+            H, psi0, tlist, c_op_mesolve, [num(2)], args
+        ).expect[0]
+
+        np.testing.assert_allclose(np.real(p_ex), np.real(p_ex_ref), atol=1e-4)
+
+    def testFloquetMasterEquation2(self):
+        """
+        Test Floquet-Markov Master Equation for a two-level system
+        subject to dissipation.
+        """
+        delta = 1.0 * 2 * np.pi
+        eps0 = 1.0 * 2 * np.pi
+        A = 0.0 * 2 * np.pi
+        omega = np.sqrt(delta ** 2 + eps0 ** 2)
+        T = (2 * np.pi) / omega
+        tlist = np.linspace(0.0, 2 * T, 101)
+        psi0 = rand_ket(2)
+        H0 = - eps0 / 2.0 * sigmaz() - delta / 2.0 * sigmax()
+        H1 = A / 2.0 * sigmax()
+        args = {'w': omega}
+        H = QobjEvo([H0, [H1, lambda t, w: np.sin(w * t)]],
+                    args=args)
+        e_ops = [num(2)]
+        gamma1 = 1
+
+        # Collapse operator for Floquet-Markov Master Equation
+        c_op_fmmesolve = sigmax()
+
+        # Collapse operator for Lindblad Master Equation
+        def spectrum(omega):
+            return (omega > 0) * omega * 0.5 * gamma1 / (2 * np.pi)
+
+        ep, vp = H0.eigenstates()
+        op0 = vp[0]*vp[0].dag()
+        op1 = vp[1]*vp[1].dag()
+
+        c_op_mesolve = _convert_c_ops(c_op_fmmesolve, spectrum, vp, ep)
+
+        # Solve the floquet-markov master equation
+        p_ex = fmmesolve(
+            H, psi0, tlist, [c_op_fmmesolve], [num(2)],
+            [spectrum], T, args=args
+        ).expect[0]
+
+        # Compare with mesolve
+        p_ex_ref = mesolve(
+            H, psi0, tlist, c_op_mesolve, [num(2)], args
+        ).expect[0]
+
+        np.testing.assert_allclose(np.real(p_ex), np.real(p_ex_ref), atol=1e-4)
+
+    @pytest.mark.parametrize("kmax", [5, 25, 100])
+    def testFloquetMasterEquation3(self, kmax):
+        """
+        Test Floquet-Markov Master Equation for a two-level system
+        subject to dissipation with internal transform of fmmesolve
+        """
+        delta = 1.0 * 2 * np.pi
+        eps0 = 1.0 * 2 * np.pi
+        A = 0.0 * 2 * np.pi
+        omega = np.sqrt(delta ** 2 + eps0 ** 2)
+        T = (2 * np.pi) / omega
+        tlist = np.linspace(0.0, 2 * T, 101)
+        psi0 = rand_ket(2)
+        H0 = - eps0 / 2.0 * sigmaz() - delta / 2.0 * sigmax()
+        H1 = A / 2.0 * sigmax()
+        args = {'w': omega}
+        H = QobjEvo([H0, [H1, lambda t, w: np.sin(w * t)]],
+                    args=args)
+        e_ops = [num(2)]
+        gamma1 = 1
+
+        # Collapse operator for Floquet-Markov Master Equation
+        c_op_fmmesolve = sigmax()
+
+        # Collapse operator for Lindblad Master Equation
+        def spectrum(omega):
+            return (omega > 0) * 0.5 * gamma1 * omega / (2 * np.pi)
+
+        ep, vp = H0.eigenstates()
+        op0 = vp[0]*vp[0].dag()
+        op1 = vp[1]*vp[1].dag()
+
+        c_op_mesolve = _convert_c_ops(c_op_fmmesolve, spectrum, vp, ep)
+
+        # Solve the floquet-markov master equation
+        floquet_basis = FloquetBasis(H, T)
+        solver = FMESolver(
+            floquet_basis, [(c_op_fmmesolve, spectrum)], kmax=kmax
+        )
+        solver.start(psi0, tlist[0])
+        p_ex = [expect(num(2), solver.step(t)) for t in tlist]
+
+        # Compare with mesolve
+        output2 = mesolve(H, psi0, tlist, c_op_mesolve, [num(2)], args)
+        p_ex_ref = output2.expect[0]
+
+        np.testing.assert_allclose(np.real(p_ex), np.real(p_ex_ref),
+                                   atol=5 * 1e-4)
+
+    def testFloquetMasterEquation_multiple_coupling(self):
+        """
+        Test Floquet-Markov Master Equation for a two-level system
+        subject to dissipation with multiple coupling operators
+        """
+        delta = 1.0 * 2 * np.pi
+        eps0 = 1.0 * 2 * np.pi
+        A = 0.0 * 2 * np.pi
+        omega = np.sqrt(delta ** 2 + eps0 ** 2)
+        T = (2 * np.pi) / omega
+        tlist = np.linspace(0.0, 2 * T, 101)
+        psi0 = rand_ket(2)
+        H0 = - eps0 / 2.0 * sigmaz() - delta / 2.0 * sigmax()
+        H1 = A / 2.0 * sigmax()
+        args = {'w': omega}
+        H = QobjEvo([H0, [H1, lambda t, w: np.sin(w * t)]],
+                    args=args)
+        e_ops = [num(2)]
+        gamma1 = 1
+
+        # Collapse operator for Floquet-Markov Master Equation
+        c_ops_fmmesolve = [sigmax(), sigmay()]
+
+        # Collapse operator for Lindblad Master Equation
+        def noise_spectrum1(omega):
+            return (omega > 0) * 0.5 * gamma1 * omega/(2*np.pi)
+
+        def noise_spectrum2(omega):
+            return (omega > 0) * 0.5 * gamma1 / (2 * np.pi)
+
+        noise_spectra = [noise_spectrum1, noise_spectrum2]
+
+        ep, vp = H0.eigenstates()
+        op0 = vp[0]*vp[0].dag()
+        op1 = vp[1]*vp[1].dag()
+
+        c_op_mesolve = []
+
+        # Convert the c_ops for fmmesolve to c_ops for mesolve
+        for c_op_fmmesolve, noise_spectrum in zip(c_ops_fmmesolve,
+                                                   noise_spectra):
+            c_op_mesolve += _convert_c_ops(
+                c_op_fmmesolve, noise_spectrum, vp, ep
+            )
+
+        # Solve the floquet-markov master equation
+        output1 = fmmesolve(
+            H, psi0, tlist, c_ops_fmmesolve, e_ops, noise_spectra, T
+        )
+        p_ex = output1.expect[0]
+        # Compare with mesolve
+        output2 = mesolve(H, psi0, tlist, c_op_mesolve, e_ops, args)
+        p_ex_ref = output2.expect[0]
+
+        np.testing.assert_allclose(np.real(p_ex), np.real(p_ex_ref), atol=1e-4)
+
+    def testFloquetRates(self):
+        """
+        Compare rate transition and frequency transitions to analytical
+        results for a driven two-level system, for different drive amplitudes.
+        """
+        # Parameters
+        wq = 5.4 * 2 * np.pi
+        wd = 6.7 * 2 * np.pi
+        delta = wq - wd
+        T = 2 * np.pi / wd
+        tlist = np.linspace(0.0, 2 * T, 101)
+        array_A = np.linspace(0.001, 3 * 2 * np.pi, 10, endpoint=False)
+        H0 = Qobj(wq / 2 * sigmaz())
+        args = {'wd': wd}
+        c_ops = sigmax()
+        gamma1 = 1
+        kmax = 5
+
+        array_ana_E0 = [-np.sqrt((delta / 2)**2 + a**2) for a in array_A]
+        array_ana_E1 = [np.sqrt((delta / 2)**2 + a**2) for a in array_A]
+        array_ana_delta = [
+            2 * np.sqrt((delta / 2)**2 + a**2)
+            for a in array_A
+        ]
+
+        def noise_spectrum(omega):
+            return (omega > 0) * 0.5 * gamma1 * omega/(2*np.pi)
+
+        idx = 0
+        for a in array_A:
+            # Hamiltonian
+            H1_p = a * sigmap()
+            H1_m = a * sigmam()
+            H = QobjEvo(
+                [H0, [H1_p, lambda t, wd: np.exp(-1j * wd * t)],
+                [H1_m, lambda t, wd: np.exp(1j * wd * t)]], args=args
+            )
+
+            floquet_basis = FloquetBasis(H, T)
+            DeltaMatrix = _floquet_delta_tensor(floquet_basis.e_quasi, kmax, T)
+            X = _floquet_X_matrices(floquet_basis, [c_ops], kmax)
+
+            # Check energies
+            deltas = np.ndarray.flatten(DeltaMatrix)
+
+            # deltas and array_ana_delta have at least 1 value in common.
+            assert (min(abs(deltas - array_ana_delta[idx])) < 1e-4)
+
+            # Check matrix elements
+            Xs = np.concatenate(
+                [X[0][k].to_array() for k in range(-kmax, kmax+1)]
+            ).flatten()
+
+            normPlus = np.sqrt(a**2 + (array_ana_E1[idx] - delta / 2)**2)
+            normMinus = np.sqrt(a**2 + (array_ana_E0[idx] - delta / 2)**2)
+
+            Xpp_p1 = (a / normPlus**2) * (array_ana_E1[idx] - delta / 2)
+            assert (min(abs(Xs - Xpp_p1)) < 1e-4)
+            Xpp_m1 = (a / normPlus**2) * (array_ana_E1[idx] - delta / 2)
+            assert (min(abs(Xs - Xpp_m1)) < 1e-4)
+            Xmm_p1 = (a / normMinus**2) * (array_ana_E0[idx] - delta / 2)
+            assert (min(abs(Xs - Xmm_p1)) < 1e-4)
+            Xmm_m1 = (a / normMinus**2) * (array_ana_E0[idx] - delta / 2)
+            assert (min(abs(Xs - Xmm_m1)) < 1e-4)
+            Xpm_p1 = (a / (normMinus * normPlus)
+                      * (array_ana_E0[idx]- delta / 2))
+            assert (min(abs(Xs - Xpm_p1)) < 1e-4)
+            Xpm_m1 = (a / (normMinus * normPlus)
+                      * (array_ana_E1[idx] - delta / 2))
+            assert (min(abs(Xs - Xpm_m1)) < 1e-4)
+            idx += 1
diff --git a/qutip/tests/solver/test_mcsolve.py b/qutip/tests/solver/test_mcsolve.py
index fdac5eb55b..f191de0d8c 100644
--- a/qutip/tests/solver/test_mcsolve.py
+++ b/qutip/tests/solver/test_mcsolve.py
@@ -385,7 +385,7 @@ def test_super_H():
     np.testing.assert_allclose(mc_expected.expect[0], mc.expect[0], atol=0.5)
 
 
-def test_McSolver_run():
+def test_MCSolver_run():
     size = 10
     a = qutip.QobjEvo([qutip.destroy(size), 'coupling'], args={'coupling':0})
     H = qutip.num(size)
@@ -405,7 +405,7 @@ def test_McSolver_run():
     assert res.num_trajectories == 1001
 
 
-def test_McSolver_stepping():
+def test_MCSolver_stepping():
     size = 10
     a = qutip.QobjEvo([qutip.destroy(size), 'coupling'], args={'coupling':0})
     H = qutip.num(size)
diff --git a/qutip/tests/solver/test_propagator.py b/qutip/tests/solver/test_propagator.py
index 6cb3ac598b..2d36d9f977 100644
--- a/qutip/tests/solver/test_propagator.py
+++ b/qutip/tests/solver/test_propagator.py
@@ -131,7 +131,7 @@ def testPropSolver(solver):
     assert (U(1.5, 0.5) - propagator(H, 1, c_ops)).norm('max') < 1e-4
 
 
-def testPropMcSolver():
+def testPropMCSolver():
     a = destroy(5)
     H = a.dag()*a
     solver = MCSolver(H, [a])