docs: def _deltasummation(term, ranges, i_range):
      """Partially execute a summation for `term` with a Kronecker Delta for one
      of the summation indices.
  
      This implements the solution to the core sub-problem in
      :func:`indexed_sum_over_kronecker`
  
      Args:
          term (QuantumExpression): term of the sum
          ranges (list): list of all summation index ranges
              (class:`IndexRangeBase` instances)
          i_range (int): list-index of element in `ranges` which should be
              eliminated
  
      Returns:
          ``(result, flag)`` where `result` is a list
          of ``(new_term, new_ranges)`` tuples and `flag` is an integer.
  
      There are three possible cases, indicated by the returned `flag`. Consider
      the following setup::
  
          >>> i, j, k = symbols('i, j, k', cls=IdxSym)
          >>> i_range = IndexOverList(i, (0, 1))
          >>> j_range = IndexOverList(j, (0, 1))
          >>> ranges = [i_range, j_range]
          >>> def A(i, j):
          ...    from sympy import IndexedBase
          ...    return OperatorSymbol(StrLabel(IndexedBase('A')[i, j]), hs=0)
  
      1. If executing the sum produces a single non-zero term, result will be
      ``[(new_term, new_ranges)]`` where `new_ranges` contains the input `ranges`
      without the eliminated range specified by `i_range`.  This should be the
      most common case for calls to:func:`_deltasummation`::
  
          >>> term = KroneckerDelta(i, j) * A(i, j)
          >>> result, flag = _deltasummation(term, [i_range, j_range], 1)
          >>> assert result == [(A(i, i), [i_range])]
          >>> assert flag == 1
  
      2. If executing the sum for the index symbol specified via `index_range`
      does not reduce the sum, the result will be the list ``[(term, ranges)]``
      with unchanged `term` and `ranges`::
  
          >>> term = KroneckerDelta(j, k) * A(i, j)
          >>> result, flag = _deltasummation(term, [i_range, j_range], 0)
          >>> assert result == [(term, [i_range, j_range])]
          >>> assert flag == 2
  
      This case also covers if there is no Kroncker delta in the term::
  
          >>> term = A(i, j)
          >>> result, flag = _deltasummation(term, [i_range, j_range], 0)
          >>> assert result == [(term, [i_range, j_range])]
          >>> assert flag == 2
  
      3. If `term` does not contain a Kronecker delta as a factor, but in a
      sum that can be expanded, the result will be a list of
      ``[(summand1, ranges), (summand2, ranges), ...]`` for the summands of that
      expansion. In this case, `:func:`_deltasummation` should be called again
      for every tuple in the list, with the same `i_range`::
  
          >>> term = (KroneckerDelta(i, j) + 1) * A(i, j)
          >>> result, flag = _deltasummation(term, [i_range, j_range], 1)
          >>> assert result == [
          ...     (A(i, j), [i_range, j_range]),
          ...     (KroneckerDelta(i,j) * A(i, j), [i_range, j_range])]
          >>> assert flag == 3
      """
      from qnet.algebra.core.abstract_quantum_algebra import QuantumExpression
      idx = ranges[i_range].index_symbol
      summands = _expand_delta(term, idx)
      if len(summands) > 1:
          return [(summand, ranges) for summand in summands], 3
      else:
          delta, expr = _extract_delta(summands[0], idx)
      if not delta:
          return [(term, ranges)], 2
      solns = sympy.solve(delta.args[0] - delta.args[1], idx)
      assert len(solns) > 0  # I can't think of an example that might cause this
      #     if len(solns) == 0:
      #         return [(term._zero, [])], 4
      if len(solns) != 1:
          return [(term, ranges)], 2
      value = solns[0]
      new_term = expr.substitute({idx: value})
      if _RESOLVE_KRONECKER_WITH_PIECEWISE:
          new_term *= ranges[i_range].piecewise_one(value)
      assert isinstance(new_term, QuantumExpression)
      return [(new_term, ranges[:i_range] + ranges[i_range+1:])], 1