# Floquet Formalism

## Introduction

Many time-dependent problems of interest are periodic. The dynamics of such systems can be solved for directly by numerical integration of the Schrödinger or Master equation, using the time-dependent Hamiltonian. But they can also be transformed into time-independent problems using the Floquet formalism. Time-independent problems can be solve much more efficiently, so such a transformation is often very desirable.

In the standard derivations of the Lindblad and Bloch-Redfield master equations the Hamiltonian describing the system under consideration is assumed to be time independent. Thus, strictly speaking, the standard forms of these master equation formalisms should not blindly be applied to system with time-dependent Hamiltonians. However, in many relevant cases, in particular for weak driving, the standard master equations still turns out to be useful for many time-dependent problems. But a more rigorous approach would be to rederive the master equation taking the time-dependent nature of the Hamiltonian into account from the start. The Floquet-Markov Master equation is one such a formalism, with important applications for strongly driven systems (see e.g., [Gri98]).

Here we give an overview of how the Floquet and Floquet-Markov formalisms can be used for solving time-dependent problems in QuTiP. To introduce the terminology and naming conventions used in QuTiP we first give a brief summary of quantum Floquet theory.

## Floquet theory for unitary evolution

The Schrödinger equation with a time-dependent Hamiltonian \(H(t)\) is

(1)\[ H(t)\Psi(t) = i\hbar\frac{\partial}{\partial t}\Psi(t),\]

where \(\Psi(t)\) is the wave function solution. Here we are interested in problems with periodic time-dependence, i.e., the Hamiltonian satisfies \(H(t) = H(t+T)\) where \(T\) is the period. According to the Floquet theorem, there exist solutions to (1) of the form

(2)\[ \Psi_\alpha(t) = \exp(-i\epsilon_\alpha t/\hbar)\Phi_\alpha(t),\]

where \(\Psi_\alpha(t)\) are the Floquet states (i.e., the set of wave function solutions to the Schrödinger equation), \(\Phi_\alpha(t)=\Phi_\alpha(t+T)\) are the periodic Floquet modes, and \(\epsilon_\alpha\) are the quasienergy levels. The quasienergy levels are constants in time, but only uniquely defined up to multiples of \(2\pi/T\) (i.e., unique value in the interval \([0, 2\pi/T]\)).

If we know the Floquet modes (for \(t \in [0,T]\)) and the quasienergies for a particular \(H(t)\), we can easily decompose any initial wavefunction \(\Psi(t=0)\) in the Floquet states and immediately obtain the solution for arbitrary \(t\)

(3)\[ \Psi(t) = \sum_\alpha c_\alpha \Psi_\alpha(t) = \sum_\alpha c_\alpha \exp(-i\epsilon_\alpha t/\hbar)\Phi_\alpha(t),\]

where the coefficients \(c_\alpha\) are determined by the initial wavefunction \(\Psi(0) = \sum_\alpha c_\alpha \Psi_\alpha(0)\).

This formalism is useful for finding \(\Psi(t)\) for a given \(H(t)\) only if we can obtain the Floquet modes \(\Phi_a(t)\) and quasienergies \(\epsilon_\alpha\) more easily than directly solving (1). By substituting (2) into the Schrödinger equation (1) we obtain an eigenvalue equation for the Floquet modes and quasienergies

(4)\[ \mathcal{H}(t)\Phi_\alpha(t) = \epsilon_\alpha\Phi_\alpha(t),\]

where \(\mathcal{H}(t) = H(t) - i\hbar\partial_t\). This eigenvalue problem could be solved analytically or numerically, but in QuTiP we use an alternative approach for numerically finding the Floquet states and quasienergies [see e.g. Creffield et al., Phys. Rev. B 67, 165301 (2003)]. Consider the propagator for the time-dependent Schrödinger equation (1), which by definition satisfies

\[U(T+t,t)\Psi(t) = \Psi(T+t).\]

Inserting the Floquet states from (2) into this expression results in

\[U(T+t,t)\exp(-i\epsilon_\alpha t/\hbar)\Phi_\alpha(t) = \exp(-i\epsilon_\alpha(T+t)/\hbar)\Phi_\alpha(T+t),\]

or, since \(\Phi_\alpha(T+t)=\Phi_\alpha(t)\),

\[U(T+t,t)\Phi_\alpha(t) = \exp(-i\epsilon_\alpha T/\hbar)\Phi_\alpha(t) = \eta_\alpha \Phi_\alpha(t),\]

which shows that the Floquet modes are eigenstates of the one-period propagator. We can therefore find the Floquet modes and quasienergies \(\epsilon_\alpha = -\hbar\arg(\eta_\alpha)/T\) by numerically calculating \(U(T+t,t)\) and diagonalizing it. In particular this method is useful to find \(\Phi_\alpha(0)\) by calculating and diagonalize \(U(T,0)\).

The Floquet modes at arbitrary time \(t\) can then be found by propagating \(\Phi_\alpha(0)\) to \(\Phi_\alpha(t)\) using the wave function propagator \(U(t,0)\Psi_\alpha(0) = \Psi_\alpha(t)\), which for the Floquet modes yields

\[U(t,0)\Phi_\alpha(0) = \exp(-i\epsilon_\alpha t/\hbar)\Phi_\alpha(t),\]

so that \(\Phi_\alpha(t) = \exp(i\epsilon_\alpha t/\hbar) U(t,0)\Phi_\alpha(0)\). Since \(\Phi_\alpha(t)\) is periodic we only need to evaluate it for \(t \in [0, T]\), and from \(\Phi_\alpha(t \in [0,T])\) we can directly evaluate \(\Phi_\alpha(t)\), \(\Psi_\alpha(t)\) and \(\Psi(t)\) for arbitrary large \(t\).

### Floquet formalism in QuTiP

QuTiP provides a family of functions to calculate the Floquet modes and quasi energies, Floquet state decomposition, etc., given a time-dependent Hamiltonian.

Consider for example the case of a strongly driven two-level atom, described by the Hamiltonian

(5)\[ H(t) = -\frac{1}{2}\Delta\sigma_x - \frac{1}{2}\epsilon_0\sigma_z + \frac{1}{2}A\sin(\omega t)\sigma_z.\]

In QuTiP we can define this Hamiltonian as follows:

    >>> delta = 0.2 * 2*np.pi
    >>> eps0 = 1.0 * 2*np.pi
    >>> A = 2.5 * 2*np.pi
    >>> omega = 1.0 * 2*np.pi
    >>> H0 = - delta/2.0 * sigmax() - eps0/2.0 * sigmaz()
    >>> H1 = A/2.0 * sigmaz()
    >>> args = {'w': omega}
    >>> H = [H0, [H1, 'sin(w * t)']]

The \(t=0\) Floquet modes corresponding to the Hamiltonian (5) can then be calculated using the `FloquetBasis` class, which encapsulates the Floquet modes and the quasienergies:

    >>> T = 2*np.pi / omega
    >>> floquet_basis = FloquetBasis(H, T, args)
    >>> f_energies = floquet_basis.e_quasi
    >>> f_energies
    array([-2.83131212,  2.83131212])
    >>> f_modes_0 = floquet_basis.mode(0)
    >>> f_modes_0
    [Quantum object: dims = [[2], [1]], shape = (2, 1), type = ket
    Qobj data =
    [[ 0.72964231+0.j      ]
     [-0.39993746+0.554682j]],
    Quantum object: dims = [[2], [1]], shape = (2, 1), type = ket
    Qobj data =
    [[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).

    >>> 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)
    >>> spsi0 = basis(2, 0)
    >>> q_energies = np.zeros((len(A_vec), 2))
    >>> H0 = delta / 2.0 * sigmaz() - eps0 / 2.0 * sigmax()
    >>> args = {'w': omega}
    >>> for idx, A in enumerate(A_vec):
    >>>   H1 = A / 2.0 * sigmax()
    >>>   H = [H0, [H1, lambda t, w: np.sin(w * t)]]
    >>>   floquet_basis = FloquetBasis(H, T, args)
    >>>   q_energies[idx,:] = floquet_basis.e_quasi
    >>> plt.figure()
    >>> plt.plot(A_vec/omega, q_energies[:,0] / delta, 'b', A_vec/omega, q_energies[:,1] / delta, 'r')
    >>> plt.xlabel(r'$A/\omega$')
    >>> plt.ylabel(r'Quasienergy / $\Delta$')
    >>> plt.title(r'Floquet quasienergies')
    >>> plt.show()

Given the Floquet modes at \(t=0\), we obtain the Floquet mode at some later time \(t\) using `FloquetBasis.mode`:

    >>> f_modes_t = floquet_basis.mode(2.5)
    >>> f_modes_t
    [Quantum object: dims = [[2], [1]], shape = (2, 1), type = ket
    Qobj data =
    [[-0.89630512-0.23191946j]
     [ 0.37793106-0.00431336j]],
    Quantum object: dims = [[2], [1]], shape = (2, 1), type = ket
    Qobj data =
    [[-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 (5) given some initial state \(\left|\psi_0\right>\). To do that, we first need to decompose the initial state in the Floquet states, using the function `FloquetBasis.to_floquet_basis`

    >>> psi0 = rand_ket(2)
    >>> f_coeff = floquet_basis.to_floquet_basis(psi0)
    >>> f_coeff
    [(-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 (5) at an arbitrary time \(t\) using the function `FloquetBasis.from_floquet_basis`:

    >>> t = 10 * np.random.rand()
    >>> 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 (5).

    import numpy as np
    from matplotlib import pyplot

    import qutip

    delta = 0.2 * 2*np.pi
    eps0  = 1.0 * 2*np.pi
    A     = 0.5 * 2*np.pi
    omega = 1.0 * 2*np.pi
    T      = (2*np.pi)/omega
    tlist  = np.linspace(0.0, 10 * T, 101)
    psi0   = qutip.basis(2, 0)

    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, w: np.sin(w * t)]]

    # 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 = floquetbasis.to_floquet_basis(psi0)

    # calculate the wavefunctions using the from the floquet modes coefficients
    p_ex = np.zeros(len(tlist))
    for n, t in enumerate(tlist):
        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
    p_ex_ref = qutip.mesolve(H, psi0, tlist, e_ops=[qutip.num(2)], args=args).expect[0]

    # plot the results
    pyplot.plot(tlist, np.real(p_ex),     'ro', tlist, 1-np.real(p_ex),     'bo')
    pyplot.plot(tlist, np.real(p_ex_ref), 'r',  tlist, 1-np.real(p_ex_ref), 'b')
    pyplot.xlabel('Time')
    pyplot.ylabel('Occupation probability')
    pyplot.legend(("Floquet $P_1$", "Floquet $P_0$", "Lindblad $P_1$", "Lindblad $P_0$"))
    pyplot.show()

### 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 times. The list of times to pre-compute modes for may be passed to `FloquetBasis` using `precompute=tlist`, and then `FloquetBasis.from_floquet_basis` and `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:

    import numpy as np
    from matplotlib import pyplot
    import qutip

    delta = 0.0  * 2*np.pi
    eps0  = 1.0 * 2*np.pi
    A     = 0.25 * 2*np.pi
    omega = 1.0 * 2*np.pi
    T      = 2*np.pi / omega
    tlist  = np.linspace(0.0, 10 * T, 101)
    psi0   = qutip.basis(2,0)

    H0 = - delta/2.0 * qutip.sigmax() - eps0/2.0 * qutip.sigmaz()
    H1 = A/2.0 * qutip.sigmax()
    args = {'w': omega}
    H = [H0, [H1, lambda t, w: np.sin(w * t)]]

    # 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 = floquetbasis.to_floquet_basis(psi0)

    # calculate the wavefunctions using the from the floquet modes coefficients
    p_ex = np.zeros(len(tlist))
    for n, t in enumerate(tlist):
        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
    p_ex_ref = qutip.mesolve(H, psi0, tlist, e_ops=[qutip.num(2)], args=args).expect[0]

    # plot the results
    pyplot.plot(tlist, np.real(p_ex),     'ro', tlist, 1-np.real(p_ex),     'bo')
    pyplot.plot(tlist, np.real(p_ex_ref), 'r',  tlist, 1-np.real(p_ex_ref), 'b')
    pyplot.xlabel('Time')
    pyplot.ylabel('Occupation probability')
    pyplot.legend(("Floquet $P_1$", "Floquet $P_0$", "Lindblad $P_1$", "Lindblad $P_0$"))
    pyplot.show()

Note that the parameters and the Hamiltonian used in this example is not the same as in the previous section, and hence the different appearance of the resulting figure.

For convenience, all the steps described above for calculating the evolution of a quantum system using the Floquet formalisms are encapsulated in the function `fsesolve`. Using this function, we could have achieved the same results as in the examples above using

    output = fsesolve(H, psi0=psi0, tlist=tlist, e_ops=[qutip.num(2)], args=args)
    p_ex = output.expect[0]

## Floquet theory for dissipative evolution

A driven system that is interacting with its environment is not necessarily well described by the standard Lindblad master equation, since its dissipation process could be time-dependent due to the driving. In such cases a rigorious approach would be to take the driving into account when deriving the master equation. This can be done in many different ways, but one way common approach is to derive the master equation in the Floquet basis. That approach results in the so-called Floquet-Markov master equation, see Grifoni et al., Physics Reports 304, 299 (1998) for details.

For a brief summary of the derivation, the important contents for the implementation in QuTiP are listed below.

The floquet mode \(\ket{\phi_\alpha(t)}\) refers to a full class of quasienergies defined by \(\epsilon_\alpha + k \Omega\) for arbitrary \(k\). Hence, the quasienenergy difference between two floquet modes is given by

\[\Delta_{\alpha \beta k} = \frac{\epsilon_\alpha - \epsilon_\beta}{\hbar} + k \Omega\]

For any coupling operator \(q\) (given by the user) the matrix elements in the floquet basis are calculated as:

\[X_{\alpha \beta k} = \frac{1}{T} \int_0^T dt \; e^{-ik \Omega t} \bra{\phi_\alpha(t)}q\ket{\phi_\beta(t)}\]

From the matrix elements and the spectral density \(J(\omega)\), the decay rate \(\gamma_{\alpha \beta k}\) is defined:

\[\gamma_{\alpha \beta k} = 2 \pi J(\Delta_{\alpha \beta k}) | X_{\alpha \beta k}|^2\]

The master equation is further simplified by the RWA, which makes the following matrix useful:

\[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})\]

The density matrix of the system then evolves according to:

\[\dot{\rho}_{\alpha \alpha}(t) = \sum_\nu (A_{\alpha \nu} \rho_{\nu \nu}(t) - A_{\nu \alpha} \rho_{\alpha \alpha} (t))\]
\[\dot{\rho}_{\alpha \beta}(t) = -\frac{1}{2} \sum_\nu (A_{\nu \alpha} + A_{\nu \beta}) \rho_{\alpha \beta}(t) \qquad \alpha \neq \beta\]

### The Floquet-Markov master equation in QuTiP

The QuTiP function `fmmesolve` implements the Floquet-Markov master equation. It calculates the dynamics of a system given its initial state, a time-dependent Hamiltonian, a list of operators through which the system couples to its environment and a list of corresponding spectral-density functions that describes the environment. In contrast to the `mesolve` and `mcsolve`, and the `fmmesolve` does characterize the environment with dissipation rates, but extract the strength of the coupling to the environment from the noise spectral-density functions and the instantaneous Hamiltonian parameters (similar to the Bloch-Redfield master equation solver `brmesolve`).

Note

Currently the `fmmesolve` can only accept a single environment coupling operator and spectral-density function.

The noise spectral-density function of the environment is implemented as a Python callback function that is passed to the solver. For example:

    gamma1 = 0.1
    def noise_spectrum(omega):
        return (omega>0) * 0.5 * gamma1 * omega/(2*pi)

The other parameters are similar to the `mesolve` and `mcsolve`, and the same format for the return value is used `Result`. The following example extends the example studied above, and uses `fmmesolve` to introduce dissipation into the calculation

    import numpy as np
    from matplotlib import pyplot
    import qutip

    delta = 0.0  * 2*np.pi
    eps0  = 1.0 * 2*np.pi
    A     = 0.25 * 2*np.pi
    omega = 1.0 * 2*np.pi
    T      = 2*np.pi / omega
    tlist  = np.linspace(0.0, 20 * T, 301)
    psi0   = qutip.basis(2,0)

    H0 = - delta/2.0 * qutip.sigmax() - eps0/2.0 * qutip.sigmaz()
    H1 = A/2.0 * qutip.sigmax()
    args = {'w': omega}
    H = [H0, [H1, lambda t, w: np.sin(w * t)]]

    # noise power spectrum
    gamma1 = 0.1
    def noise_spectrum(omega):
        return (omega>0) * 0.5 * gamma1 * omega/(2*np.pi)

    # solve the floquet-markov master equation
    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_coeff_t = output.floquet_states[idx]
        psi_t = output.floquet_basis.from_floquet_basis(f_coeff_t, t)
        # Alternatively
        psi_t = output.states[idx]
        p_ex[idx] = qutip.expect(qutip.num(2), psi_t)

    # For reference: calculate the same thing with mesolve
    output = qutip.mesolve(
        H, psi0, tlist, [np.sqrt(gamma1) * qutip.sigmax()],
        e_ops=[qutip.num(2)], args=args
    )
    p_ex_ref = output.expect[0]

    # plot the results
    pyplot.plot(tlist, np.real(p_ex), 'r--', tlist, 1-np.real(p_ex), 'b--')
    pyplot.plot(tlist, np.real(p_ex_ref), 'r', tlist, 1-np.real(p_ex_ref), 'b')
    pyplot.xlabel('Time')
    pyplot.ylabel('Occupation probability')
    pyplot.legend(("Floquet $P_1$", "Floquet $P_0$", "Lindblad $P_1$", "Lindblad $P_0$"))
    pyplot.show()

Finally, `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()], e_ops=[num(2)],
                       spectra_cb=[noise_spectrum], T=T, args=args)
    p_ex = output.expect[0]