\section{Supplemental Materials}

In this Supplement we provide a discussion of the detectability of gravitational waves
from stellar-mass binary black holes (BBHs) by the Laser Interferometer Space
Antenna (LISA), an alternative to the multiparameter test presented in the paper  and
the accuracy of Fisher matrix inversion.

\paragraph*{Archival searches for stellar-mass BBHs in the LISA data:}
Ground-based observatories such as the Cosmic Explorer (CE) \cite{2019arXiv190704833R} 
and Einstein Telescope (ET) \cite{Punturo:2010zz}, have the best sensitivity to 
stellar-mass BBHs of 5-100 $M_\odot$ and can detect them up to redshifts of $z\sim10$ 
and beyond. A small fraction of such BBHs that are close enough will also be observable 
by LISA \cite{LISA2017}. The observability of a source depends on the false alarm rate at a 
given signal-to-noise-ratio (SNR) and the number of trials needed to dig out the signal. 

As an example, there is only a chance of 0.13\% (i.e., a p-value of 0.0013)
that a single draw from a Gaussian distribution with zero mean and unit variance 
yields a number larger than 3.  On the other hand, multiple draws from the same 
distribution increases the p-value for getting a number larger than 3. 
Likewise, if it is necessary to carry out a 
blind search for GW signals without any knowledge of their parameters then one ought to 
employ a large number of templates which makes it computationally expensive 
to dig out weaker signals 
from a noisy background.  Indeed, Ref.~\cite{Moore:2019pke} argued that as many
as $10^{40}$ templates would be needed to dig out a stellar-mass BBH signal from 
LISA data. This would require a matched filter SNR $\geq14$ for a p-value of 
$10^{-3}$. 

A third-generation (3G) network of CE and ET can aid in searching for stellar-mass BBH signals in LISA 
data as the former will detect them with extremely high fidelity and can therefore provide 
a very tight prior on the source's extrinsic (sky position, polarization, orientation 
of the orbit, and luminosity distance) as well as intrinsic (masses and spins of the 
companion black holes) parameters. Therefore, an archival search for signals can greatly 
reduce the number of templates/trials required and hence enhance their detectability 
in the LISA data.

The analyses in Ref.\,\cite{Moore:2019pke} assumed for the LISA archival search 
that the trigger time, the epoch when the signal's amplitude reaches its maximum
value in the detector, is known precisely, while the errors in the intrinsic and other 
extrinsic parameters will be smaller in 3G detectors than those measured by
LIGO \cite{Discovery} by a factor of 10. This contracted volume of the search space
decreases the number of templates by a factor of $10^{29}$ compared to 
a blind search \cite{Moore:2019pke}.  Yet the result is that LISA would need 
$\sim10^{11}$ templates to identify a stellar-mass BBH signal, or an SNR 
threshold of $\sim$9 for a p-value of $10^{-3}.$ Other studies \cite{Wong:2018uwb} 
have shown that using an alternate method one can detect GW150914-like 
source even with an SNR$\sim$7 for the same p-value.  We have computed the 
number of templates for the events considered in this work and find that LISA 
archival searches require far fewer templates owing to the greatly contracted
search volume thanks to the high-fidelity of 3G observations.

\begin{figure*}
\includegraphics[width=0.98\textwidth]{hist-supp.pdf}
\caption{The figure shows histograms of several parameters for the best 95\% of the events
that are visible both in 3G ground-based detectors and LISA (with an SNR $\ge 4$):
(i) the SNR, (ii) the error in the chirpmass, (iii) the error in 
the sky position, all determined by the 3G network, and (iv) the number of templates 
required to search for a stellar-mass BBH in the LISA data.}
\label{fig:num templates}
\end{figure*}

In one year, the ground-based 3G network of CE and ET would find ${\sim}200$ events 
that would have an SNR of 4 or more in LISA. For each such event we estimated 
the number of templates necessary for an archival search in LISA data. To this end, we 
computed the Fisher information matrix and its inverse---the variance-covariance matrix---to 
obtain the expected error in the measurement of various extrinsic and intrinsic
parameters using a 3G network of two CE detectors, one each in the USA and 
Australia, and one ET in Europe. Fig.\,\ref{fig:num templates} plots a subset of our
results: the distribution of the SNR, the error in the measurement of the chirpmass, 
and the uncertainty in the sky position (the left three panels). In making these plots
we have rejected the worst 5\% of the outliers in our simulation.

Due to the large SNR in the 3G network the parameters of the events are constrained 
very tightly. We assume that the events are contained within the 2-sigma region of the 
parameter uncertainties. This means that 90\% of the sources 
will be resolved by the 3G network to within $\sim 60$ arcmin$^2,$ which is better than 
the angular resolution of LISA for these sources; at an SNR of $\rho=10$, 
$\Delta\Omega_{\rm LISA} \sim [1.22 \lambda/(D\rho)]^2\sim 1800\,\rm arcmin^2,$ 
where $D=2\,\rm AU$ is LISA's baseline and $\lambda\sim3\times 10^7\,\rm m$ is the wavelength of GWs
corresponding to a frequency of 10 mHz.  
Moreover, the angles describing the orientation of the detector are also precisely 
determined by the 3G network---errors in the inclination and polarization angles are measured 
to within 1.0$^\circ$ and 2.3$^\circ$, respectively, for 90\% of the events.  Consequently, it is not 
necessary to include the angular parameters in the search nor the trigger time, which will be known
to better than 15 $\mu$s for 90\% of the events. For the Gaussian distribution of spins considered
in this paper it is not necessary to include spins either, leaving just the two masses over 
which a search should be carried out.

The number of templates required for an archival search for stellar-mass BBHs in the 
\emph{LISA band} for the ${\sim}200$ events in our population is shown in the right most 
panel of Fig.\,\ref{fig:num templates}. These numbers are computed using a minimal match 
of 0.95 (or allowing for a loss of less than 5\% in the SNR) \cite{Owen:1998dk}
and assuming that the true event lies in the 2-sigma region of the uncertainty in the masses
determined by the 3G network.  A vast majority (90\%) of the events require $< 15,000$ templates. 
One can employ singular value decomposition to find the number of \emph{independent} templates 
\cite{Cannon:2011vi}, which is typically 1 to 2 orders-of-magnitude smaller than the 
number of templates found at this minimal match, or ${\sim}150$-$1500$ for most of the events. 
Thus, an SNR-4 event in LISA will have a p-value $10^{-2}$ or smaller. 


\paragraph*{Multiparameter tests from the higher PN side:}
One may consider an interesting variant of the multiparameter tests of GR where more than one PN parameter is treated as independent, starting from the highest order that is currently known, which is 3.5PN. This may be thought of as tests of alternatives to GR where up to a particular PN order, the predictions of both GR and its alternative match but beyond that they deviate. This may be naturally motivated from an effective field theoretic perspective where the deviations may appear when the binary dynamics proceeds beyond a certain scale of velocity or field strength~\cite{Endlich:2017tqa,Sennett:2019bpc}. 

\begin{figure}[hb]
\includegraphics[scale=0.35]{right_SNR4_new_with_GB.pdf}
	\caption{Multiparameter tests using multiband observations with LISA and CE. 
         Shown are combined $1\sigma$ bounds on various PN coefficients starting 
         from 3.5PN through 0PN in the inspiral phase of the signal while measuring 
         many of them together at a time. 
	  Different types of markers symbolize how 
         many PN deformation parameters were constrained simultaneously. 
        For example, `$\CIRCLE$'  represents `one PN deformation parameter at a time’, 
         `$\vee$' represents `two PN deformation parameters at a time’, and so on.
         The filled diamonds and pentagons are bounds obtained with CE and LISA,
         respectively, on the last four and five PN coefficients from their 
         respective golden binaries.  The total masses of the CE 
         and LISA golden binaries are $200 \,M_{\odot}$ 
         and $6.6\times10^{5}\,M_{\odot}$, both binaries are 1 Gpc away and have 
         component spins $\chi_1=0.6, \chi_2=0.5$.
         The figure represents results for the BBH population having Gaussian spin 
	 distribution (we get similar estimates for a uniform spin distribution).
	 }\label{fig:boundsright}
\end{figure}

Figure~\ref{fig:boundsright} shows the 
results for multiparameter tests starting from 3.5PN through 0PN order.
 For instance a three parameter test, would correspond to the case where only the last three PN deformation parameters (3.5PN, 3PN {\it log} and 3PN) are simultaneously estimated and so on.
The most significant result here is for the 7-parameter test for
which it is found that  the last seven PN parameters  can be bounded with
$\lesssim 10\%$ using the population we simulated earlier. 
This means, if we assume that the leading Newtonian 
coefficient is not modified, the simultaneous constraints on the remaining
seven are of the order of a few percent. Since a modification to the 
Newtonian phasing  would be at odds with the extremely stringent
bounds on them from binary pulsar observations~\cite{Yunes:2010qb}, this is the most general test we wish to carry out from this perspective.
Fig.~\ref{fig:boundsright} also shows the best bounds that CE  and LISA alone could yield from their golden binaries.
This demonstrates that multibanding of GW signals 
is the only way to put meaningful constraints on multiple PN deformation parameters with higher accuracy.


\paragraph*{Inversion Accuracy of the Fisher matrix:}
Due to large degeneracies among the PN deformation parameters, there are high chances that the
Fisher matrices (obtained from different events in the simulated population) corresponding to different multiparameter tests will be rendered ill-conditioned. 
Such Fisher matrices when inverted might lead to unreliable bounds. We therefore impose an inversion accuracy criterion on the Fisher matrices 
for all multiparameter tests. This criterion is defined to be $|\Gamma\cdot\Sigma -{\rm I}|\leq {\cal
O}(10^{-3})$, where $\Gamma$, $\Sigma$, and ${\rm I}$ are the multiband multiparameter Fisher matrix, the
corresponding variance-covariance matrix and the identity matrix, respectively. Any Fisher matrix obtained from an event corresponding a 
particular multi-parameter test, that does not meet this criterion is dropped from our analysis and is not used to obtain the combined multiband multiparameter bounds. 



\bibliographystyle{apsrev}
