The phase interaction models for disperse multiphase flows are crucial for modeling various scientific and engineering applications. In the existing theory of disperse multiphase flows, phase interactions are represented as forces between the phases. It is well known that models using only force-based interactions are ill-posed, resulting in non-hyperbolic governing equations[1]. The multiphase flow community generally agrees that the issue of ill-posedness arises from the absence of correct physics, which hinders numerical stability and convergence in simulations. In this work, using ensemble averaging and nearest particle statistics, the fluid-particle phase interaction is decomposed into a particle-mean-field force (representing the conventional drag force) and the divergence of the particle-fluid-particle (PFP) stress[2,3]. To explore the physics of PFP stress and develop a practical mathematical model, a series of particle-resolved simulations have been conducted using high-order computational fluid dynamics-discrete element method coupled by the immersed boundary method. In our work[1], we found that for a fluid-particle system with high particle-fluid density ratios, the PFP stress is attractive along the flow direction and repulsive in the direction perpendicular to the flow. This serves as the macroscopic representation of the well-known drafting-kissing-tumbling (DKT) phenomena observed in experiments[4], which is still absent in existing particle-laden flow theories. To quantify the essential physics, a preliminary PFP stress model, expressed in terms of particle volume fraction and Reynolds number, is proposed based on statistical mechanics and particle-resolved direct numerical simulations. Recent studies[5] show that our PFP stress model not only aligns the numerical results of particle-laden flows with complex experimental observations but also makes the governing equations of particle-laden flows hyperbolic — a significant breakthrough for the theory of disperse multiphase flows.