Explore Trotterisation of the super-Linbladian to simulate noisy dynamics

[TysonRayJones]

Dynamical evolution under a time-independent Hamiltonian and bath-interaction solves the master equation as

|rho(t)> = exp(L t) |rho(0)>

where L is a superoperator form of the Linbladian (containing the Hamiltonian and jump operators), and |rho> vectorises the density matrix.

See Sec 3.2.2 "The Markovian quantum master equation" of Breuer and Petruccione's 'The Theory of Open Quantum Systems' helpfully linked by Truman Yu Ng!

The term exp(L t) can be Trotterised into a series of tractable super-operators which QuEST can individually simulate. This permits noisy simulation without integrating the differential master equation.

Explore whether this is robust and suitable for QuEST! One could imagine interfaces like

applyUnitaryTimeEvolution(
    Qureg qureg, PauliStrSum hamiltonian, qreal time, int order, int reps);

applyImaginaryTimeEvolution(
    Qureg qureg, PauliStrSum hamiltonian, qreal tau, int order, int reps);

applyNoisyTimeEvolution(
    Qureg qureg, PauliStrSum hamiltonian, CompMatr* jumps, int numJumps, qreal time, int order, int reps);

where the first two wrap (for convenience) the existing apply(NonUnitary)TrotterizedPauliStrSumGadget() and the final uses an extension of the internal Trotter functions to "mix-in" non-Pauli-gadgets (the matrix exponential of the jumps operators must be found). We could alternatively restrict to only jump operators expressed in the Pauli basis (limited but still interesting), express L entirely as a non-Hermitian weighted sum of Pauli strings, and trotterise exclusively with the existing functions.

[TysonRayJones]

Trivial when jump operators are diagonal or Pauli strings but otherwise requires exponentiating arbitrary, dense matrices. Potentially implement one or more schemes here after determining reasonable assumptions of jump operators
https://arxiv.org/abs/2404.12789

[TysonRayJones]

A functional demo accepting only PauliStr jump operators is here 🎉 But it remains prudent to also support projectors (which also admit optimised treatment), diagonal matrices (optimisable too!) and dense matrices (requiring some nuisance).

A function which simultaneously accepts all of these input formats (PauliStr, list of projector outcomes, DiagMatr, CompMatr and possibly more) is no small order. A possible solution is to make a new abstracting JumpOp struct which stores each format in separate fields. E.g.

JumpOp createJumpOpFromPauliStr(PauliStr);
JumpOp createJumpOpFromCompMatr(CompMatr);
JumpOp createJumpOpFromDiagMatr(DiagMatr);
JumpOp createJumpOpFromMultiQubitProjector(outcomes, numOutcomes);

// needs more thought
JumpOp createJumpOpFromPauliStrSum(PauliStrSum);
JumpOp createJumpOpFromOffDiagonalMultiQubitProjector(braOutcomes, ketOutcomes, numOutcomes);

These can be passed independent of their target qubits (so that the jumps can be re-targeted) via

int numJumps = 3;
qreal jumpRates[] = {.4, .1, 0.8};
JumpOp jumpOps[] = {j1, j2, j3};
int jumpTargs[][] = {
    {1, 2},    // must match dimension of j1
    {0, 1, 2}, // ... j2
    {2}        // ... j3
};

applyTrotterizedNoisyTimeEvolution(
    qureg, hamil, 
    jumpRates, jumpOps, jumpTargs, numJumps, 
    time, order, reps
);

Then applyTrotterizedNoisyTimeEvolution() can handle each JumpOp optimally for its underlying representation. This design also enables a few sneaky optimisations;

• when given a CompMatr, JumpOp will pre-allocate the needed matrices for storing exp(a J* (x) J), exp(a dag(J) . J) and exp(a transp(J) . conj(J)). Though the latter two are conjugates of one another, pre-computing and storing both may be worthwhile to avoid superfluous conj calls in every Trotter repetition.
• when applyTrotterizedNoisyTimeEvolution() is invoked, CompMatr-type JumpOp must evaluate the aforementioned matrix exponentials using the passed correspondingjumpRate and time. These "witnessed variables" can be saved to the JumpOp to avoid recalculation of the aforementioned matrices when subsequent invocation uses the same damping rates and time.

[TysonRayJones]

A functional demo accepting jump operators as PauliStrSum is here 🎉

[TysonRayJones]

That function (applyTrotterizedNoisyTimeEvolution) was committed to devel (6eec1cf), though an extension to the additional jump operator forms described above remains possible.