def series_expand(
            self, param: Symbol, about, order: int) -> tuple:
        r"""Expand the expression as a truncated power series in a
        scalar parameter.

        When expanding an expr for a parameter $x$ about the point $x_0$ up to
        order $N$, the resulting coefficients $(c_1, \dots, c_N)$ fulfill

        .. math::

            \text{expr} = \sum_{n=0}^{N} c_n (x - x_0)^n + O(N+1)

        Args:
            param: Expansion parameter $x$
            about (Scalar): Point $x_0$ about which to expand
            order: Maximum order $N$ of expansion (>= 0)

        Returns:
            tuple of length ``order + 1``, where the entries are the
            expansion coefficients, $(c_0, \dots, c_N)$.

        Note:
            The expansion coefficients are
            "type-stable", in that they share a common base class with the
            original expression. In particular, this applies to "zero"
            coefficients::

                >>> expr = KetSymbol("Psi", hs=0)
                >>> t = sympy.symbols("t")
                >>> assert expr.series_expand(t, 0, 1) == (expr, ZeroKet)
        """
        expansion = self._series_expand(param, about, order)
        # _series_expand is generally not "type-stable", so we continue to
        # ensure the type-stability
        res = []
        for v in expansion:
            if v == 0 or v.is_zero:
                v = self._zero
            elif v == 1:
                v = self._one
            assert isinstance(v, self._base_cls)
            res.append(v)
        return tuple(res)
r"""Expand the expression as a truncated power series in a
        scalar parameter.

        When expanding an expr for a parameter $x$ about the point $x_0$ up to
        order $N$, the resulting coefficients $(c_1, \dots, c_N)$ fulfill

        .. math::

            \text{expr} = \sum_{n=0}^{N} c_n (x - x_0)^n + O(N+1)

        Args:
            param: Expansion parameter $x$
            about (Scalar): Point $x_0$ about which to expand
            order: Maximum order $N$ of expansion (>= 0)

        Returns:
            tuple of length ``order + 1``, where the entries are the
            expansion coefficients, $(c_0, \dots, c_N)$.

        Note:
            The expansion coefficients are
            "type-stable", in that they share a common base class with the
            original expression. In particular, this applies to "zero"
            coefficients::

                >>> expr = KetSymbol("Psi", hs=0)
                >>> t = sympy.symbols("t")
                >>> assert expr.series_expand(t, 0, 1) == (expr, ZeroKet)