Molecular composition in the interstellar medium (ISM) is rich: so far, more than 300 molecular species have been identified through their spectroscopic signatures. If we also include the polycyclic aromatic hydrocarbons (PAHs), which are estimated to contain up to 10% of all galactic carbon, the potential molecular variety grows even larger. However, the ISM is still a harsh place for molecules: cosmic radiation from various sources provides high-energy photons that can provide single and multiple ionization of the molecular species. For instance, if a PAH molecule will absorb an extreme UV or soft X-ray photon, this can lead up to triple ionization with a production of a stable trication and/or to fragmentation. Therefore, molecular mono-and polycations can exist in the ISM, serving as intermediates in different chemical reactions and generating new molecular species. The potential role of multiply charged ions in the production of complex molecules can be crucial since the rates of ion-neutral reactions grow linearly with the increase of the ionic charge. Over the years, many experimental results were obtained in investigating molecular polycations. For instance, dications were detected in space, in the tails of comets, trications of PAHs were produced with high-energy photons, by above-threshold ionization (ATI), and also by interaction with highly-charged xenon ions. The quest for molecular tetracations started with the high-power ATI by femtosecond intense laser pulses, which led to the laboratory detection of tetracations of PAHs, substituted benzene, and even doubly-iodosubstituted acetylene. Intact tetracations of PAHs were also found in high-energy collisions of neutral PAH molecules with various ions. However, to our knowledge, no existing model provides a general description of the stability of polycations. In the case of fragmentation, the availability of many possible pathways of fragmentation possibilities with various spin-charge states provides a very challenging chemical space to sample. Theoretical calculations may be done, for instance, using molecular dynamics-based approaches, reproducing some of the polycationic properties, such as the kinetic energy release in fragmentation. Semi-quantitative models are known to be useful for describing trends of different processes. In particular, a model for polycation stability based on the Coulomb interaction potential was proposed in Ref. , and for fullerenes, a Hückel-method based model was introduced in Ref. .