When an extreme shock wave releases from the free surface of the material, some high speed particulate matters will be ejected from the material body and enter into the background gas. This induced multiphase mixing phenomenon is known as the ejecta mixing. Ejecta mixing is one of the most important problems in the scope of inner explosive compression engineering, and it is also a frontier research subject of the impact dynamics, multiphase fluid dynamics, computational mathematics, etc. The properties of ejecta mixing have been investigated experimently and analytically for many years. However, the results of numerical simuliation are very rare. At present, the ejecta mixing study mainly focuses on the gas particle one-way coupling, that is, the interests of existing works are in the characteristics of the ejected particulate matter transport in the gas. In fact, after a large number of particles entering into the gas, the gas and the particles will interact with each other. So it is necessary to consider the feedback of particles to the gas. In this paper, the theoretical modeling of gas particle two-way coupling, the discrete algorithm of the mathematical model and the particle phase feedback effects on the gas shock wave propagation are investigated in the framework of Lagrangian coordinates. In order to obtain the details of ejecta movement, the particle trajectory model is chosen as the basic model, and then the governing equations including interactions between gas phase and particle phase are derived. For giving the specific calculation formula, the physical meanings of the coupled interaction source terms in the Lagrangian framework are analyzed and a stable numerical scheme is given based on the staggered strategy. We also devise two different computing models of ejecta mixing, the planar and the column configurations, and then the numerical simulations are carried out. The phenomenon of gas shock speed acceleration caused by particle feedback is found, and the distributions of the physical quantities, such as density, velocity, specific internal energy, pressure, in the gas area are changed. Especially for the convergent configuration, the feedback effects will be amplified further by the geometrical shrinkage, which may have a significant influence on the performance of the inner explosion compression, owing to the obvious uniformity variation of the gas flow field and the gas shock rebound in advance. The mathematical model, the numerical method and the new physical findings in this paper, will provide an important theoretical support for the in-depth understanding of the ejecta mixing and also for the solving of the corresponding engineering application problems.