docs: def factor_for_trace(ls: HilbertSpace, op: Operator) -> Operator:
      r'''Given a :class:`.LocalSpace` `ls` to take the partial trace over and an
      operator `op`, factor the trace such that operators acting on disjoint
      degrees of freedom are pulled out of the trace. If the operator acts
      trivially on ls the trace yields only a pre-factor equal to the dimension
      of ls. If there are :class:`LocalSigma` operators among a product, the
      trace's cyclical property is used to move to sandwich the full product by
      :class:`LocalSigma` operators:
  
      .. math::
  
          {\rm Tr} A \sigma_{jk} B = {\rm Tr} \sigma_{jk} B A \sigma_{jj}
  
      Args:
          ls: Degree of Freedom to trace over
          op: Operator to take the trace of
  
      Returns:
          The (partial) trace over the operator's spc-degrees of freedom
      '''
      if op.space == ls:
          if isinstance(op, OperatorTimes):
              pull_out = [o for o in op.operands if o.space is TrivialSpace]
              rest = [o for o in op.operands if o.space is not TrivialSpace]
              if pull_out:
                  return (OperatorTimes.create(*pull_out) *
                          OperatorTrace.create(OperatorTimes.create(*rest),
                                               over_space=ls))
          raise CannotSimplify()
      if ls & op.space == TrivialSpace:
          return ls.dimension * op
      if ls < op.space and isinstance(op, OperatorTimes):
          pull_out = [o for o in op.operands if (o.space & ls) == TrivialSpace]
  
          rest = [o for o in op.operands if (o.space & ls) != TrivialSpace]
          if (not isinstance(rest[0], LocalSigma) or
              not isinstance(rest[-1], LocalSigma)):
              for j, r in enumerate(rest):
                  if isinstance(r, LocalSigma):
                      m = r.j
                      rest = (
                          rest[j:] + rest[:j] +
                          [LocalSigma.create(m, m, hs=ls), ])
                      break
          if not rest:
              rest = [IdentityOperator]
          if len(pull_out):
              return (OperatorTimes.create(*pull_out) *
                      OperatorTrace.create(OperatorTimes.create(*rest),
                                           over_space=ls))
      raise CannotSimplify()