docs: def move_drive_to_H(slh, which=None, expand_simplify=True):
      r'''Move coherent drives from the Lindblad operators to the Hamiltonian.
  
      For the given SLH model, move inhomogeneities in the Lindblad operators (resulting
      from the presence of a coherent drive, see :class:`CoherentDriveCC`) to the
      Hamiltonian.
  
      This exploits the invariance of the Lindblad master equation under the
      transformation  (cf. Breuer and Pettrucione, Ch 3.2.1)
  
      .. math::
          :nowrap:
  
          \begin{align}
              \Op{L}_i &\longrightarrow \Op{L}_i' = \Op{L}_i - \alpha_i  \\
              \Op{H}   &\longrightarrow
              \Op{H}' = \Op{H} + \frac{1}{2i} \sum_j
                      (\alpha_j \Op{L}_j^{\dagger} - \alpha_j^* \Op{L}_j)
          \end{align}
  
      In the context of SLH, this transformation is achieved by feeding `slh` into
  
      .. math::
          \SLH(\identity, -\mat{\alpha}, 0)
  
      where $\mat{\alpha}$ has the elements $\alpha_i$.
  
      Parameters
      ----------
      slh : SLH
          SLH model to transform. If `slh` does not contain any inhomogeneities, it is
          invariant under the transformation.
  
      which : sequence or None
          Sequence of circuit dimensions to apply the transform to. If None, all
          dimensions are transformed.
  
      expand_simplify : bool
          if True, expand and simplify the new SLH object before returning. This has no
          effect if `slh` does not contain any inhomogeneities.
  
      Returns
      -------
      new_slh : SLH
          Transformed SLH model.
  
      '''
      if which is None:
          which = []
      scalarcs = []
      for jj, L in enumerate(slh.Ls):
          if not which or jj in which:
              scalarcs.append(-get_coeffs(L.expand())[IdentityOperator])
          else:
              scalarcs.append(0)
  
      if np.all(np.array(scalarcs) == 0):
          return slh
      new_slh = SLH(identity_matrix(slh.cdim), scalarcs, 0) << slh
      if expand_simplify:
          return new_slh.expand().simplify_scalar()
      return new_slh