docs: def detunings_combinations(pairs):
      r"""Return all combinations of detunings.
  
      >>> Ne = 6
      >>> Nl = 2
      >>> omega_level = [0.0, 100.0, 100.0, 200.0, 200.0, 300.0]
      >>> xi = np.zeros((Nl, Ne, Ne))
      >>> coup = [[(1, 0), (2, 0)], [(3, 0), (4, 0), (5, 0)]]
      >>> for l in range(Nl):
      ...     for pair in coup[l]:
      ...         xi[l, pair[0], pair[1]] = 1.0
      ...         xi[l, pair[1], pair[0]] = 1.0
  
      >>> aux = define_simplification(omega_level, xi, Nl)
      >>> u, invu, omega_levelu, Neu, xiu = aux
      >>> pairs = detunings_indices(Neu, Nl, xiu)
      >>> detunings_combinations(pairs)
      [[(1, 0), (2, 0)], [(1, 0), (3, 0)]]
  
      """
      def iter(pairs, combs, l):
          combs_n = []
          for i in range(len(combs)):
              for j in range(len(pairs[l])):
                  combs_n += [combs[i] + [pairs[l][j]]]
          return combs_n
  
      Nl = len(pairs)
      combs = [[pairs[0][k]] for k in range(len(pairs[0]))]
      for l in range(1, Nl):
          combs = iter(pairs, combs, 1)
  
      return combs