The excess of muons in ultra-high energy cosmic rays (UHECR) in comparison with the results of Monte-Carlo simulations that grows with the primary particle energy is known as the muon puzzle. Since the LHC data became available, many improvements have been done to the models of hadronic interactions used in extensive air showers (EAS) simulation tools. Yet the noticeable deviation from experimental data is observed for both so-called post-LHC models: EPOS LHC and QGSJET-II-04. The focus of theoretical works that try to solve this puzzle has shifted towards the process of hadronization, i. e. when partons produced in high-energy collisions combine into final state hadrons. Some attempts were made to improve the description of particle production like core-corona effect in EPOS. But for the complete solution, further changes are required.
This study presents the detailed overview of models of hadronization used in several modern high-energy hadron interaction models: PYTHIA 8.3, EPOS LHC and QGSJET-II-04. The connection to the EAS simulations and muon puzzle is considered. Different possibilities for further improvements of the models are discussed.