We investigate the intrinsic and observational properties of $z\gtrsim 6$
galaxies hosting coalescing massive black holes (MBHs) that gives rise to
gravitational waves (GWs) detectable with the Laser Interferometer Space
Antenna (LISA). We adopt a zoom-in cosmological hydrodynamical simulation of
galaxy formation and black hole (BH) co-evolution, zoomed-in on a $M_h \sim
10^{12}~\rm M_{\odot}$ dark matter halo at z = 6, which hosts a fast accreting
super-massive black hole (SMBH) and a star-forming galaxy. Following the SMBH
formation backward in time, we identify the merging events that concurred to
its formation and we pick up the ones that are detectable with LISA. Among
these LISA detectable events (LDEs), we select those that, based on their
intrinsic properties are expected to be bright in one or more electromagnetic
(EM) bands. We post-process these events with dust radiative transfer
calculations to make predictions about their spectral energy distributions and
continuum maps in the JWST to ALMA wavelength range. We compare the spectra
arising from galaxies hosting the merging MBHs with those arising from AGN
powered by single accreting BHs. We find that it will be impossible to identify
an LDE from the continuum SEDs because of the absence of specific imprints from
the merging MBHs. We also compute the profile of the H$_{\rm \alpha}$ line
arising from LDEs, considering the contribution from their star-forming regions
and the accreting MBHs. We find that the presence of two accreting MBHs would
be difficult to infer even if both MBHs accrete at super-Eddington rates. We
conclude that the combined detection of GW and EM signals from $z\gtrsim 6$
MBHs is challenging not only because of the poor sky-localization provided by
LISA, but also because the loudest GW emitters are not massive enough to leave
significant signatures in the emission lines arising from the broad line
region.
