We employ theoretically “exact” and numerically “accurate” Beyond Born–Oppenheimer (BBO) treatment to construct diabatic potential energy surfaces (PESs) of the benzene radical cation (C6H6+) for the first time and explore the workability of the time-dependent discrete variable representation (TDDVR) method for carrying out dynamical calculations to evaluate the photoelectron (PE) spectra of its neutral analog. Ab initio adiabatic PESs and nonadiabatic coupling terms are computed over a series of pairwise normal modes, which exhibit rich nonadiabatic interactions starting from Jahn–Teller interactions and accidental conical intersections/seams to pseudo Jahn–Teller couplings. Once the electronic structure calculation is completed on the low-lying five doublet electronic states (X̃2E1g, B̃2E2g, and C̃2A2u) of the cationic species, diabatization is carried out employing the adiabatic-to-diabatic transformation (ADT) equations for the five-state sub-Hilbert space to compute highly accurate ADT angles, and thereby, single-valued, smooth, symmetric, and continuous diabatic PESs and couplings are constructed. Subsequently, such surface matrices are used to perform multi-state multi-mode nuclear dynamics for simulating PE spectra of benzene. Our theoretical findings clearly depict that the spectra for X̃2E1g and B̃2E2g−C̃2A2u states obtained from BBO treatment and TDDVR dynamics exhibit reasonably good agreement with the experimental results as well as with the findings of other theoretical approaches.