code: def outer_product(vec0: QubitVector, vec1: QubitVector) -> QubitVector:
      
      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)
docs: Direct product of qubit vectors
  
      The tensor ranks must match and qubits must be disjoint.