docs: def detunings_indices(Neu, Nl, xiu):
      r"""Get the indices of the transitions of all fields.
  
      They are returned in the form
      [[(i1, j1), (i2, j2)], ...,[(i1, j1)]].
      that is, one list of pairs of indices for each field.
  
      >>> 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
      >>> detunings_indices(Neu, Nl, xiu)
      [[(1, 0)], [(2, 0), (3, 0)]]
  
      """
      pairs = []
      for l in range(Nl):
          ind = []
          for iu in range(Neu):
              for ju in range(iu):
                  if xiu[l, iu, ju] == 1:
                      ind += [(iu, ju)]
          pairs += [ind]
      return pairs