docs: def symbolic_master_equation(self, rho=None):
          """Compute the symbolic Liouvillian acting on a state rho
  
          If no rho is given, an OperatorSymbol is created in its place.
          This correspnds to the RHS of the master equation
          in which an average is taken over the external noise degrees of
          freedom.
  
          Args:
              rho (Operator): A symbolic density matrix operator
  
          Returns:
              Operator: The RHS of the master equation.
  
          """
          L, H = self.L, self.H
          if rho is None:
              rho = OperatorSymbol('rho', hs=self.space)
          return (-I * (H * rho - rho * H) +
                  sum(Lk * rho * adjoint(Lk) -
                      (adjoint(Lk) * Lk * rho + rho * adjoint(Lk) * Lk) / 2
                      for Lk in L.matrix.ravel()))