Adding Floquet and some useful infrastructure

[Gavin-Rockwood]

Problem Description

It would be nice to have built in tools for Floquet analysis similar to qutip along with some of the functionality needed to make this easy (such as a propagator function and an ensemble ability for sesolve for easy parallelization, useful for calculating propagators).

Proposed Solution

There should be a structure Floquet Basis which contains the floquet quasi-energies and a function to return the floquet modes at a time t. This would require the implementation of a function to compute the U(ti, tf) propagator, a function that would benefit from an extension to sesolve that allows it to take in lists of wavefunctions and generate an ensemble problem for the sciml ODE solver.

Alternate Solutions

No response

Additional Context

I have had my own implementation that I am currently working on incorporating into a fork of QuantumToolbox (here. It consists of

1. The introduction of a version of sesolve that dispatches on vectors of ket Quantum objects (here). This generates an EnsembleTimeEvolutionProblem instead of a TimeEvolutionProblem (though it does contain one) and uses this to generate an EnsembleProblem for DifferentialEquations. This allows for easy multi-threading or distributed solving for a set of wavefunctions. This is used to speed up the calculation of propagators.
   • This function also has the ability to perform this ensemble solve over a list of parameters. By setting iterate_params=true, the function will treat the params kwarg as a collection of parameter sets and solve all of those as well. In the end, this function returns a $m\times n$ matrix of TimeEvolutionSol objects where the rows are the initial wavefunctions and the columns are the different sets of parameters.
2. The introduction of a function propagator which takes in a Hamiltonian "H" and a times ti, tf and passes this Hamiltonian along with the H(ti) eigenstates to sesolve to compute the propagator $\sum_i\ket{i,t=t_i}\bra{i,t=t_f}$ (here).
   - There is also a structure Propagator (constructed by get_propagator) which can be evaluated as Propagator(tf) to return the $t=t_f$ propagator (via the function .eval which Propagator contains). This is really just a convenience thing as it allows one to just create some propagator object "U" and then keep mindlessly reusing it.
3. Finally, in here I have implemented the some basic floquet tools. These include:
   1. A structure FloquetBasis, initialized by get_floquet_basis which includes
      • T: The floquet period
      • e_quasi: The list of floquet quasienergies
      • modes: A function which when evaluated as modes(t) returns the floquet modes at time "t" by using the function propagate_floquet_modes. This is done by calculating the propagator $U(t_i=0, t_f%T)$ and acting it on the $t=0$ floquet modes. If modes is evaluated at $t=0$, it skips this entirely and just returns the $t=0$ floquet modes. This function, can take in kwargs that are passed through to the propagator calculator and onto sesolve.
   2. A function `floquet_sweep' which takes in a function "H(params)" and a list params to sample and calculates the a FloquetBasis for each parameter in the list and passes this list of states and spectra to track_stater which sorts the states across all the parameters. This is useful for doing stuff like searching for resonances by calculating the floquet spectrum for a range of drive frequencies and searching for avoided crossings. There are probably better algorithms to do this but this one that I used to do the floquet analysis for the drive dynamics in this paper along with doing some transmon ionization analysis.

[ytdHuang]

Hi @Gavin-Rockwood,

Thank you for the idea of implementing Floquet-related methods.

As you probably know, one of the main purposes of our package is to align with qutip (Python). That is, the behavior (including function names, return values, arguments, keyword arguments, and even their order) should be as similar as possible to here. This similarity may allow users to copy-and-paste code from Python to Julia with minimal modifications and have it execute correctly.

I think the idea you proposed can be divided into the following steps:

1. Add Propagator structure, propagator function, and corresponding functionalities (which is also implemented in qutip [click here])
2. Define structure FloquetBasis and corresponding utility functions
3. Add fsesolve (solving Schrodinger equation using the Floquet formalism)
4. Add floquet_tensor and fmmesolve (Floquet-Markov master equation)

Since this is your first time contributing to this package, I think it is better that you start from a smaller pull request, which gets you understand more of the structure of our package before doing many lines of changes.

For a smaller pull request but also relate to your topic, I suggest you can implement the Propagator stuff first.

[Gavin-Rockwood]

Hi @ytdHuang,
Yeah, implementing the propagator stuff is absolutely something I can work while making sure it is inline with qutip. I know that in qutip, the propagator is calculated by time evolving an identity operator, but sesolve is not set up to do this. Would it still be useful to include the extension to sesolve that generates ensemble problems? (And I guess same for mesolve for propagators with c_ops.)

[ytdHuang]

If I understand correctly, both Propagator and propagator in qutip do deal with sesolve, where the returned propagator are just unitary operators. The propagator calculated by evolving an identity operator is just an alternative way for other solvers. See here for more details.

[ytdHuang]

For the ensemble problems, you just want to multi-thread over different initial conditions (states) and maybe different parameters of Hamiltonian right ?

I think this can be done totally from user's side with Julia's built-in multi-threading feature.

If we really want to do it from our side, then I expect it will be a huge change of our time evolution solvers in order to remain type-stability.

Therefore, I think we should put this idea on hold, and consider it in the future roadmap. But maybe @albertomercurio has different opinion ?

[Gavin-Rockwood]

I'm still learning the details of type stability, but I do have (what i believe, if I am reading @code_warntype correctly) a type stable version of extension of sesolve for doing ensemble problems on lists of wave functions. It is able to infer that it will be returning a vector of TimeEvolutionSol objects. Type stability issues start arising when I add multiple parameters into the mix as that is depending on a kwarg variable.

Doing multi-threading or distributed for parameters could totally be done from the users side (though the way I have it set up, its easy to define a EnsembleTimeEvolutionProblem for both initial states and parameters and call sesolve on this as well). I think my main interest with this idea is how to efficiently calculate propagators, especially for larger dynamical systems. If I do it by time evolving individual, especially if the Hamiltonian is sparse, then this would speed it up without much thought on the user side. I'm not sure however how these ensemble problems will work with CUDA. On mac using Metal, there are limitations on the types of ODEs that can be done using ensemble.

[ytdHuang]

@albertomercurio is on vacation now. We can wait for his reply about this topic.

Maybe you can first try to implement the Propagator stuff ?

[albertomercurio]

Hello everyone. Sorry I was on vacation (they are quit long in Italy).

I think the idea is really nice, and a lot of things seem to be implemented, so I would first coordinate on the steps to do. Maybe splitting into several pull requests.

Regarding the ensemble ODE for sesolve, I agree that this can speed up the generation of the propagator. However, evolving the identity can also be very fast if the Hamiltonian is on GPU. The best solution would be to have both methods available, also to align with QuTiP, where they support also operators as input for sesolve.

Beware that a kind of propagator is already implemented in the Arnoldo-Lindblad eigensolver eigsolve_al:

• https://github.com/qutip/QuantumToolbox.jl/blob/main/src/qobj/eigsolve.jl#L109