Limits of vacuum-template subtraction for LISA massive black hole binary sources in realistic environments

We investigate the impact of gravitational wave (GW) dephasing due to gas accretion on the subtraction of massive black hole (MBH) binary signals over 4 yr of LISA data in the context of the global-fit. Based on state of the art predictions for the population of merging MBHs, we show that imperfect subtraction with vacuum waveform templates leaves a GW residual with an SNR of $3.2^{+5.4}_{-1.9}\times \sqrt{f_{\rm Edd} \langle \dot n \rangle/(20\, {\rm yr}^{-1})}$, where $f_{\rm Edd}$ is the typical Eddington ratio and $\langle \dot n \rangle$ the mean merger rate of LISA MBH binaries. We characterize the dependence of the residual on key population hyper-parameters and provide a simple fitting function. {Finally, we consider the distinguishability of the residual as a stochastic signal by comparing with Bayesian power-law sensitivity curves, while discussing several additional detection and mitigation strategies. Overall, our analysis indicates that the residual is unlikely to be confidently distinguishable from instrumental noise while nevertheless being likely to bias the inference of other signals.