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

    Parameters
    ----------
    self : Builder or Approach
        The system given as Builder or Approach object.
    b_,bp_ : int
        Labels of the many-body states.

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