  In the next decade, the Laser Interferometer Space Antenna (LISA) will detect
the coalescence of massive black hole binaries (MBHBs) in the range $[10^4,
10^8] \, \rm M_{\odot}$, up to $z\sim10$. Their gravitational wave (GW) signal
is expected to be accompanied by an electromagnetic counterpart (EMcp),
generated by the gas accreting on the binary or on the remnant BH. In this
work, we present the number and characteristics (such as redshift and mass
distribution, apparent magnitudes or fluxes) of EMcps detectable jointly by
LISA and some representative EM telescopes. We combine state-of-the-art
astrophysical models for the galaxies formation and evolution to build the
MBHBs catalogues, with Bayesian tools to estimate the binary sky position
uncertainty from the GW signal. Exploiting additional information from the
astrophysical models, such as the amount of accreted gas and the BH spins, we
evaluate the expected EM emission in the soft X-ray, optical and radio bands.
Overall, we predict between 7 and 21 EMcps in 4 yrs of joint observations by
LISA and the considered EM facilities, depending on the astrophysical model. We
also explore the impact of the hydrogen and dust obscuration of the optical and
X-ray emissions, as well as of the collimation of the radio emission: these
effects reduce the number to EMcps to 2 or 3, depending on the astrophysical
model, again in 4 yrs of observations. Most of the EMcps are characterised by
faint EM emission, challenging the observational capabilities of future
telescopes. Finally, we also find that systems with multi-modal sky position
posterior distributions represent only a minority of cases and do not affect
significantly the number of EMcps.
