Combined time evolution of quantum, classical, stochastic, etc. systems

[david-pl]

A more general interface to pass quantum systems in combination with others to an ODE solver (for example DifferentialEquations.jl) would be nice. This is currently a bit tedious since it requires ugly modifications of modules (see e.g. a semi-classical system where the classical differential equations are stochastic here). It's probably better to this after #92.

[ChrisRackauckas]

Can you describe what exactly you're looking for?

[david-pl]

We're looking to have arbitrary combinations of quantum systems and classical systems as we already do with the semiclassical module, but now also allowing for stochastic terms on either or both (quantum and classical) parts.

So basically, we can add a method for schroedinger_dynamic and master_dynamic, where when a function fstoch is specified, it spawns an SDEProblem rather than an ODEProblem. This can identically be done for semiclassical systems. We might have to make an additional distinction there, though, to specify whether the noise applies to the quantum or classical part (something like fstoch_quantum, fstoch_classical).

Stochastic Schrödinger and master equations #152 describing e.g. homodyne detection are then a special case, where a specific form of noise affects the quantum part.

[ChrisRackauckas]

That sounds good. A very similar integrate! front-end needs to then be added to handle stochastic, or it needs to be extended a bit. All that's really required is a quick wrapper to recast! to states before calling the user-derivatives and then adding the callback which does saving via fout. But it's so similar that it may just want to type-check and add a few things in order to re-use most of the code. Let me know if you want help here.

[david-pl]

@ChrisRackauckas I have started implementing stochastic equations in the stochastic barnch. So far only the Schrödinger equation is implemented. Before continuing I wanted to ask if you could take a glance at the code and tell me whether I got things right.

To clarify what I did so far:

Basically I want to implement an equation of the form dpsi = H*psi*dt + Hs*psi*dW. The function stochastic.schroedinger_dynamic takes two functions fdeterm and fstoch, which return the deterministic and stochastic operators H and Hs, respectively. Each of these functions is then mapped to vector form using recast! and passed as the functions f and g to an SDEProblem.

So far so good. Now, what I am not too sure about is the way I implemented a stochastic Schrödinger equation that allows for multiple Wiener increments, so simply dpsi = H*psi*dt + sum(Hs_i*psi*dW_i). I allowed for fstoch to return a Tuple of operators. Then for each operator in the tuple the product Hs_i*psi is calculated. This is then returned as an 1xn array of the resulting Kets, where n is the number of operators in the tuple. That array is the output of the g function of the SDEProblem. Hence, the solver should calculate the sum of the Kets in the return array, where each Ket is multiplied by a distinct Wiener increment. Is that correct?

If you have time, could you also take a look at the integrate_stoch function in timeevolution_base.jl? Is SRA1() a good default algorithm choice for these types of problems?

Of course any comments regarding performance, type stability, clean syntax, etc. are also very welcome. Thanks!

[david-pl]

Two things to note:

• The stochastic Schrödinger equation in its standard form has to be interpreted as Stratonovich SDE. Hence, the choice of SRA1 as default algorithm is wrong. I had some decent success in reproducing known physical examples with RKMil(interpretation=:Stratonovich).
• I did not set noise_rate_prototype when there is more than one stochastic operator given. The calculation still runs through without any problems, but I have yet to find an easily comparable physical example to check whether the results are actually correct. My guess would be that as it stands now, actually the sum over the stochastic terms sum(Hs_i*psi*dW_i) is not calculated in the way I want. It might be more something in the way of sum(Hs_i*psi)*dW.

[david-pl]

This is implemented in #204