We implement a recursive Green function method to extract the Fock space (FS) propagator and associated self-energy across the many-body localization (MBL) transition, for one-dimensional interacting fermions in a random onsite potential. We show that the typical value of the imaginary part of the local FS self-energy, \Delta_t, related to the decay rate of an initially localized state, acts as a probabilistic order parameter for the thermal to MBL phase transition; and can be used to characterize critical properties of the transition as well as the multifractal nature of MBL states as a function of disorder strength W. In particular, we show that a fractal dimension D_s extracted from \Delta_t jumps discontinuously across the transition, from D_s<1 in the MBL phase to D_s= 1 in the thermal phase. Moreover, \Delta_t follows an asymmetrical finite-size scaling form across the thermal-MBL transition, where a non-ergodic volume in the thermal phase diverges with a Kosterlitz-Thouless like essential singularity at the critical point W_c, and controls the continuous vanishing of \Delta_t as W_c is approached. In contrast, a correlation length ({\xi}) extracted from \Delta_t exhibits a power-law divergence on approaching W_c from the MBL phase.
