def get_phi0(sys, b_, bp_):
    """
    Get the reduced density matrix element corresponding to
    many-body states b and bp.

    Parameters
    ----------
    sys : Approach, Approach2vN, or Builder
        The system given as Approach, Approach2vN, or Builder object.
    b\\_,bp\\_ : int
        Labels of the many-body states.

    Returns
    --------
    phi0bbp : complex
        A matrix element of the reduced density matrix (complex number).
    """
    b = sys.si.states_order[b_]
    bp = sys.si.states_order[bp_]
    bcharge = sum(sys.si.get_state(b))
    bpcharge = sum(sys.si.get_state(bp))
    phi0bbp = 0.0
    if sys.funcp.kerntype == 'Pauli':
        if b == bp:
            ind = sys.si.get_ind_dm0(b, b, bcharge, maptype=1)
            phi0bbp = sys.phi0[ind]
    elif bcharge == bpcharge:
        ind = sys.si.get_ind_dm0(b, bp, bcharge, maptype=1)
        conj = sys.si.get_ind_dm0(b, bp, bcharge, maptype=3)
        if ind != -1:
            if type(sys.si).__name__ == 'StateIndexingDMc':
                phi0bbp = sys.phi0[ind]
            else:
                ndm0, npauli = sys.si.ndm0, sys.si.npauli
                phi0bbp = sys.phi0[ind] + 1.0j * sys.phi0[ndm0 - npauli + ind
                    ] * (+1 if conj else -1) * (0 if ind < npauli else 1)
    return phi0bbp
