This thesis explores theoretical astrophysics of massive, sub-parsec black hole binaries (BHBs) in a gaseous environment with an additional focus on identifying potential multimessenger signatures. First, the context of massive and supermassive BHBs in cosmological structure formation is given and observational techniques both for electromagnetic (EM) and gravitational waves (GW) are outlined. Then, the numerical methods of (i) Newtonian smoothed particle hydrodynamics (SPH) and (ii) thermal, optically thick general relativistic radiation hydrodynamics (GR-RHD) are detailed. In the case of the GR-RHD formulation, in which stiff source terms arise, a dedicated numerical treatment in form of an implicit explicit (IMEX) Runge-Kutta time integrator is shown. This IMEX scheme is implemented, validated and tested in two GR codes with the aim that, once publicly released, it will be useful for the community. Also to this end, the physical limitations of the current framework's applicability are explored finding that the inclusion of radiation gives order of magnitude improvements over evolutions without radiation even in situations where parts of the underlying assumptions are violated. The remaining shortcomings are discussed and strategies for future improvement are listed. The astrophysical accretion scenarios employed are transonic spherical accretion and sub- and supersonic Bondi-Hoyle-Lyttleton accretion onto a single black hole of galactic size. The secular dynamics of massive BHBs in self gravitating circumbinary accretion discs is studied at the evolutionary stage where the binary has already excavated a cavity inside the disc. Starting with prograde binary-disc systems, the torques coming from the disc-gravity and from accretion are dissected into components in both space and frequency. Alternative treatments of the thermodynamics inside the cavity are tested to check the sensitivity of the numerical results towards such idealised treatments. It is explained why there is a limiting eccentricity ecrit, which a BHB attains during disc migration, whose existence is independent of the thermodynamical treatment of the cavity gas and independent of accretion torques. Ongoing work on BHBs in retrograde, self gravitating discs is presented hinting at the existence of yet another limiting eccentricity: an instability is observed where the BHB instead of growing maximally eccentric to e 1, starts tilting with respect to its disc and the angular momentum vectors of the BHB and the disc subsequently align. Focussing on the observational implications of the prograde disc scenario, the residual eccentricity at band-entrance of LISA is calculated and found to be non-negligible. Furthermore, the expected differences in eccentricity populations produced by stellar and gaseous environments is estimated finding that eccentricity can be used together with expected spin populations as a discriminant between the two environmental histories of BHBs that are detected via GWs only. For Pulsar Timing sources, analogous calculations are presented and moreover, co-incident electromagnetic counterparts in the X-ray are identified as realistically detectable with current and upcoming X-ray facilities. Specifically, periodic variability stemming from the highly eccentric BHB could be observed by MAXI or eROSITA, whereas a significant detection of double FeKa lines requires an Athena like X-ray spectrometer. These results have impact on low frequency GW data analysis, on searches for supermassive BHBs, on the development of more realistic numerical treatments of relativistic accretion flows and some of the results on secular evolution of binaries in discs could possibly be applied to stellar and planetary migration astrophysics.