Direct product of qubit vectors

    The tensor ranks must match and qubits must be disjoint.

def outer_product(vec0: QubitVector, vec1: QubitVector) -> QubitVector:
    """Direct product of qubit vectors

    The tensor ranks must match and qubits must be disjoint.
    """
    R = vec0.rank
    R1 = vec1.rank

    N0 = vec0.qubit_nb
    N1 = vec1.qubit_nb

    if R != R1:
        raise ValueError('Incompatibly vectors. Rank must match')

    if not set(vec0.qubits).isdisjoint(vec1.qubits):
        raise ValueError('Overlapping qubits')

    qubits: Qubits = tuple(vec0.qubits) + tuple(vec1.qubits)

    tensor = bk.outer(vec0.tensor, vec1.tensor)

    # Interleave (super)-operator axes
    # R = 1  perm = (0, 1)
    # R = 2  perm = (0, 2, 1, 3)
    # R = 4  perm = (0, 4, 1, 5, 2, 6, 3, 7)
    tensor = bk.reshape(tensor, ([2**N0] * R) + ([2**N1] * R))
    perm = [idx for ij in zip(range(0, R), range(R, 2*R)) for idx in ij]
    tensor = bk.transpose(tensor, perm)

    return QubitVector(tensor, qubits)