We investigate the collective dynamics of a network comprising two populations of globally coupled phase oscillators with intrinsic frequency heterogeneity and varying fractions of pairwise and higher-order interactions. Our results show that, with homogeneous phase lag parameters, increasing the fraction of higher-order interactions and coupling strength leads to more complex dynamics, including distinct monostable and bistable chimera regions. Considering the heterogeneity of the phase lag parameter between pairwise and higher-order interactions, our study reveals that increasing the fraction of higher-order interactions leads to the emergence of various bistable and multistable regions while destabilizing monostable chimera regions, especially at small coupling strengths. Conversely, increasing the coupling strength has minimal impact on the system's dynamics for small fractions of higher-order interactions, whereas a larger fraction of higher-order interactions uncovers additional bistable and multistable regions. We derive low-dimensional reduced equations from the N-dimensional discrete system using the Ott-Antonsen ansatz and obtain bifurcation curves using XPPAUT software. Additionally, we deduce stability conditions for both synchronized and desynchronized states, which align precisely with the numerical results.