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. 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 responsible for the formation of the unique chemical fingerprint of P-rich stars. We performed a detailed 1D LTE 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 lies on the neutron-capture elements (Z>30), in particular on the elements between Sr and Ba, as well as Pb, as they provide valuable constraints to nucleosynthesis calculations. In past works, we showed that the heavy elements 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-process. 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 kind of stars. Constraining upper limits could also be estimated for CdI, InI, and SnI. 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. Regarding the nucleosynthetic origin of the pattern, no global solution explaining all four stars could be found. The model that produces the least number of discrepancies in three of the four stars is a combination of s- and i-process, but the current lack of extensive multi-dimensional hydrodynamic simulations to follow the occurrence of the i-process in different types of stars makes this scenario highly uncertain.