Analysis of post-HF wavefunctions
There are! First off, you're correct that Hartree-Fock and DFT produce effective one-particle states $\phi_i$ which are the molecular orbitals (MOs) in the case of HF or the Kohn-Sham (KS) orbitals in the case of DFT (I'll focus on Hartree-Fock since I am more familiar with it than DFT). The MOs are defined via expansion coefficients ($C_{\mu i}$) in the basis set $\{\chi_\mu\}$ such that $\phi_i = \sum_\mu C_{\mu i} \chi_{\mu}$. The more physically meaningful quantity one works with is the one-particle reduced density matrix (1PRDM) defined as $\Gamma_{pq} = \sum_{\mu \nu} C^*_{\mu p} P_{\mu \nu} C_{\nu q}$ where $P_{\mu \nu} = \sum_i C_{\mu i}^*C_{\nu i}$. Diagonalization of the 1PRDM produces what are known as natural atomic orbitals (NAOs) $\phi^{\rm{NAO}}_i = \sum_p a_{p i} \phi_p$ where the coefficients $a_{p i}$ are defined by the eigenvalue relation $\Gamma_{p q} a_{q i} = \rho_i a_{p i}$ where $\rho_i$ is the occupation number corresponding to th $i^{\rm{th}}$ NAO. Since $\Gamma$ is Hermitian, the NAOs are orthonormal and form a perfectly suitable basis to perform any quantum chemistry calculation, however, their value lies in the fact that NAOs actually show you where the electron density is so it's a very useful visualization tool. Once you have NAOs, all of the analysis you want is available. There are things called natural bond orbitals (NBOs) which rotate NAOs in such a way that they maximize the electron density between 2 atoms, thus simulating the notion of a chemical bond. Here, you find the orbitals describing your $\sigma$- and $\pi$-bonds. The point I'm trying to make is that what's important for this kind of analysis is having $\Gamma_{pq}$.
So the answer to your question depends on whether or not you calculate $\Gamma_{pq}$ for a variety of post-HF methods, such as configuration interaction (CI), coupled-cluster (CC), and others. If we denote the Hartree-Fock determinant as $|\Phi\rangle$, the 1PRDM from before is simply $\Gamma_{pq} = \langle \Phi | a^+_p a_q | \Phi \rangle$. The wavefunction originating from post-HF calculations is $|\Psi_\mu\rangle = (1 + C_\mu)|\Phi\rangle$ in the case of CI and $|\Psi_\mu\rangle = R_\mu e^T |\Phi\rangle$ in the case of equation-of-motion (EOM) CC. The analogous expression of the 1PRDM for the $\mu^{\rm{th}}$ excited state ($\mu = 0$ is the ground state) is
$$\Gamma^{\mu}_{pq} = \langle \Psi_\mu | a^+_p a_q | \Psi_\mu \rangle$$
You can calculate this object directly for the various types of wavefunction that come out of correlated electronic structure calculations. It's not light work for sure, but if you know your way around many-body algebra (e.g. Slater's Rules, diagrammatics, Wick's Theorem), you can calculate this object and diagonalize it to find correlated bonding orbitals just as you would in Hartree-Fock. Also note that in the case of CC, the bra state $\langle \Psi_\mu | \neq [|\Psi_\mu\rangle]^+$ and instead we use the biorthogonal left-CC parameterization $\langle \Psi_\mu | = \langle \Phi | L_\mu e^{-T}$ where $L_\mu$ is a linear de-excitation operator complementary to $R_\mu$. These and other technicalities are of course well-documented in the CC literature. The key observation here is that we never explicitly build a correlated $|\Psi_\mu\rangle$ defined via Slater determinants directly. In the case of CC wavefunctions, you'd blow up your computer since $e^T$ is an object as large as the full CI wavefunction for any level of truncation. So the way things are done is you obtain $|\Psi_\mu\rangle$ indirectly though a set of CI coefficients or CC cluster amplitudes and calculate the 1PRDM, again, completely defined through this sequence of finite coefficients. It may still be a single-particle orbital, but it has many-body information embedded within it.
That being said, these calculations are available in most electronic structure programs.