Support for direct products

[pmenczel]

Problem Description

Currently, there is no great support for working with direct product spaces in qutip. For that reason, both the HEOM solver and qutip-qoc contain code manipulating data layer matrices directly.

Most commonly, direct product spaces arise because of coupled differential equations. For example, we might have two states $|\psi\rangle$ and $|\psi'\rangle$ satisfying

$$ \partial_t \begin{pmatrix} |\psi\rangle \\ |\psi'\rangle \end{pmatrix} = -i \begin{pmatrix} H & 0 \\ H' & H \end{pmatrix} \begin{pmatrix} |\psi\rangle \\ |\psi'\rangle \end{pmatrix} $$

then the object $( |\psi\rangle, |\psi'\rangle )$ is a direct product.

Some applications:

• HEOM is a coupled differential equation for the state $\rho$ of an open system, and a number of auxiliary density operators $\rho_{\mathbf n}$.
• The goat algorithm in qutip-qoc works by integrating coupled Schroedinger equations for two ore more unitary operators.
• In full counting statistics, we can find the first n derivatives of a moment-generating function by solving coupled Lindblad equations for n+1 operators.
• A direct product state could potentially be used instead of a tensor product state in cases where we know that two systems are only coupled classically and there will never be any coherences between them.

Proposed Solution

I think that we would need the following:

• A way to construct objects such as $( |\psi\rangle, |\psi'\rangle )$ and the matrix in the equation above from their constituents.
• A way to describe the dimensions of these objects.
• A way to extract the constituents from these objects at the end of the calculation.
• The differential equations mentioned in the applications above should be solvable with sesolve and/or mesolve.
• Note that in the first three applications, we actually have direct products of operator spaces (which might make things complicated?)

Alternate Solutions

If all the involved Hilbert spaces have the same dimension, then $\mathcal H \times \cdots \times \mathcal H \simeq \mathbb C^n \otimes \mathcal H$. In other words, a "big matrix" can be obtained by tensoring a "small matrix" with another appropriate matrix, for example

$$ \begin{pmatrix} H & 0 \\ H' & H \end{pmatrix} = \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix} \otimes H + \begin{pmatrix} 0 & 0 \\ 1 & 0 \end{pmatrix} \otimes H' . $$

This can be written in qutip already now, of course. However, proper support for this kind of object would be nicer. Also, it is not clear to me how one would apply this tensoring-approach, for example, to the coupled differential equations for unitary operators in goat.

Additional Context

No response

[ajgpitch]

Just to add that we also have operator evolution in the form

$$ \partial_t \begin{pmatrix} U \\ U' \end{pmatrix} = -i \begin{pmatrix} H & 0 \\ H' & H \end{pmatrix} \begin{pmatrix} U \\ U' \end{pmatrix}, $$

and similar with $U$ switched for a map and $H$ for some Lindblad superoperator, in goat, which seems to work fine in terms of the multiplication during the ODE update step. So, as suggested, it would seem that mainly we just want a nice, efficient method for constructing / extracting these objects.