docs: def lindblad_operator(A, rho):
      r"""This function returns the action of a Lindblad operator A on a density\
   matrix rho. This is defined as :
  
      .. math::
          \mathcal{L}(A, \rho) = A \rho A^\dagger -
          (A^\dagger A \rho + \rho A^\dagger A)/2.
  
      >>> rho=define_density_matrix(3)
      >>> lindblad_operator( ketbra(1,2,3) ,rho )
      Matrix([
      [   rho22, -rho12/2,        0],
      [-rho21/2,   -rho22, -rho23/2],
      [       0, -rho32/2,        0]])
  
      """
      a, b = sigma_operator_indices(A)
      # print(111, a, b)
      if a is not None and b is not None:
          Ne = A.shape[0]
          L = zeros(Ne, Ne)
          L[a, a] += rho[b, b]
  
          for j in range(Ne):
              L[b, j] += -rho[b, j]/2
  
          for i in range(Ne):
              L[i, b] += -rho[i, b]/2
  
          return L
      else:
          return A*rho*A.adjoint() - (A.adjoint()*A*rho + rho*A.adjoint()*A)/2