docs: def prepare_adiabatic_limit(slh, k=None):
      """Prepare the adiabatic elimination on an SLH object
  
      Args:
          slh: The SLH object to take the limit for
          k: The scaling parameter $k \rightarrow \infty$. The default is a
              positive symbol 'k'
  
      Returns:
          tuple: The objects ``Y, A, B, F, G, N``
          necessary to compute the limiting system.
      """
      if k is None:
          k = symbols('k', positive=True)
      Ld = slh.L.dag()
      LdL = (Ld * slh.L)[0, 0]
      K = (-LdL / 2 + I * slh.H).expand().simplify_scalar()
      N = slh.S.dag()
      B, A, Y = K.series_expand(k, 0, 2)
      G, F = Ld.series_expand(k, 0, 1)
  
      return Y, A, B, F, G, N