Context. The atmospheres of phosphorus-rich (P-rich) stars have been shown to contain between 10 and 100 times more P than our Sun. Given its crucial role as an essential element for life, it is especially necessary to uncover the origin of P-rich stars to gain insights into the still unknown nucleosynthetic formation pathways of P in our Galaxy. Aims. Our objective is to obtain the extensive chemical abundance inventory of four P-rich stars, covering a large range of heavy (Z > 30) elements. This characterization will serve as a milestone for the nuclear astrophysics community to uncover the processes that form the unique chemical fingerprint of P-rich stars. Methods. We performed a detailed 1D local thermodynamic equilibrium abundance analysis on the optical UVES spectra of four P-rich stars. The abundance measurements, complemented with upper-limit estimates, included 48 light and heavy elements. Our focus lay on the neutron-capture elements (Z > 30), in particular, on the elements between Sr and Ba, as well as on Pb, as they provide valuable constraints to nucleosynthesis calculations. In past works, we showed that the heavy-element observations from the first P-rich stars are not compatible with either classical s-process or r-process abundance patterns. In this work, we compare the obtained abundances with three different nucleosynthetic scenarios: a single i-process, a double i-process, and a combination of s- and i-processes. Results. We have performed the most extensive abundance analysis of P-rich stars to date, including the elements between Sr and Ba, such as Ag, which are rarely measured in any type of stars. We also estimated constraining upper limits for Cd I, In I, and Sn I. We found overabundances with respect to solar in the s-process peak elements, accompanied by an extremely high Ba abundance and slight enhancements in some elements between Rb and Sn. No global solution explaining all four stars could be found for the nucleosynthetic origin of the pattern. The model that produces the least number of discrepancies in three of the four stars is a combination of s- and i-processes, but the current lack of extensive multidimensional hydrodynamic simulations to follow the occurrence of the i-process in different types of stars makes this scenario highly uncertain.