docs: def symbolic_heisenberg_eom(
              self, X=None, noises=None, expand_simplify=True):
          """Compute the symbolic Heisenberg equations of motion of a system
          operator X.  If no X is given, an OperatorSymbol is created in its
          place.  If no noises are given, this correspnds to the
          ensemble-averaged Heisenberg equation of motion.
  
          Args:
              X (Operator): A system operator
              noises (Operator): A vector of noise inputs
  
          Returns:
              Operator: The RHS of the Heisenberg equations of motion of X.
          """
          L, H = self.L, self.H
  
          if X is None:
              X = OperatorSymbol('X', hs=(L.space | H.space))
  
          summands = [I * (H * X - X * H), ]
          for Lk in L.matrix.ravel():
              summands.append(adjoint(Lk) * X * Lk)
              summands.append(-(adjoint(Lk) * Lk * X + X * adjoint(Lk) * Lk) / 2)
  
          if noises is not None:
              if not isinstance(noises, Matrix):
                  noises = Matrix(noises)
              LambdaT = (noises.adjoint().transpose() * noises.transpose()).transpose()
              assert noises.shape == L.shape
              S = self.S
              summands.append((adjoint(noises) * S.adjoint() * (X * L - L * X))
                              .expand()[0, 0])
              summand = (((L.adjoint() * X - X * L.adjoint()) * S * noises)
                         .expand()[0, 0])
              summands.append(summand)
              if len(S.space & X.space):
                  comm = (S.adjoint() * X * S - X)
                  summands.append((comm * LambdaT).expand().trace())
  
          ret = OperatorPlus.create(*summands)
          if expand_simplify:
              ret = ret.expand().simplify_scalar()
          return ret