Rethinking search signatures: Hadronic decays of GeV-scale feebly coupled particles

Experimental searches for GeV-scale feebly interacting particles commonly target the simplest decay signatures with two charged particles, including hadrons. Hadronic modes often dominate the inclusive decay rate, and their simulation relies extensively on $\text{textsc}Pythia$. Its fragmentation model, calibrated mainly on LEP data, is unreliable in this mass range: it populates symmetry-forbidden states and misidentifies the dominant allowed channels. We construct a simple hadronization model tuned to electromagnetic scattering data and constrained by conservation laws and available exclusive decay calculations, while retaining a common prescription applicable to different particle models. Considering the SHiP experiment as an example, we show that our model shifts the dominant search signature after event selection from two-particle to multiparticle decays, including mixed charged-neutral states. We provide the model in a form that can be readily integrated into $\text{textsc}Pythia$-based experimental simulation frameworks.