This work reports the structure, intermolecular forces, electronic/optical properties, and stability in solution of complexes formed between polycyclic aromatic hydrocarbons (PAH) and phosphorene nanoflakes by density functional theory modeling. PAH molecules reach a strong affinity with phosphorene by forming well-ordered domains, whose interaction strength decreases 13–21% compared to carbonaceous surfaces, e.g., graphene. The adsorption energies are linearly related to the NH:NC ratio of PAHs, where NH and NC are the numbers of H and C atoms; consequently, the cohesive energy of phosphorene-graphene heterostructures is estimated in 44 meV/atom. Energy decomposition (ALMO-EDA) and electron-density-based analyses support the major role of electrostatics driving forces in the interaction mechanism, which is balanced with dispersion effects for larger PAHs. In addition, phosphorene-PAH complexes display outstanding stability in solution under polar/non-polar solvents, which is due to the high polarity of the complexes and strong overcompensation of destabilizing solvation energies with stabilizing electrostatic effects. Moreover, PAHs behave as n-dopants for phosphorene, inducing small bandgap opening and weak effects on the photophysical fingerprint of phosphorene. Nevertheless, strong electron acceptor/donor and larger PAHs (NH:NC < 0.5) lead to major effects on the bandgap control, acting as active sites for orbital-controlled interactions. These findings serve as a framework for further investigation of phosphorene-based materials for remediation of PAH pollutants in water treatment technologies and uses of PAHs for phosphorene surface passivation or bandgap engineering for sensing.