Question: For gravitational waves then, what actually determines the angular uncertainty of the source direction? Does it turn out to be baseline-limited (e.g. λ/D) or instrument-limited, or limited by the very short-lived and chirped nature of the signal? Or is it something else, perhaps in the modeling and reconstruction of the event itself?
In short, it's both.
Gravitational wave interferometers are all-sky detectors. Generally, a parameter-estimation pipeline is done in a Bayesian framework to coherently analyze data from all detectors in the network. The Bayesian approach maps out the posterior prob-
ability density function of the signal parameters given the data and model providing the statistical measurement uncertainty on parameter estimates, see here. However, a full Bayesian analysis is very computationally costly, so some approximate techniques have been developed. Fast and accurate ways of computing the sky location are desirable because we want to perform follow-up searches with electromagnetic telescopes to look for possible electromagnetic counterparts which is a very "time sensitive" issue.
There are two approximate methods for determining the source location on the sky of the Earth: timing triangulation and Fisher information matrix inversion. These two methods can generally produce different uncertainties, see e.g. here.
In timing triangulation (see section II.B), which requires a network of detectors (i.e., more than one), the sky location of a source can be constructed through the difference in time of arrival of GW radiation at the different detector sites. For the first detected event GW150914, Virgo was not online yet but it's sky location was still found in a ring on the sky from the 6.3 ms time delay between the two LIGO detectors' observations. The network localization accuracy is a function of the timing uncertainty in each detector and the pairwise separation vectors of the detectors in the network. A two network detector only provides partial localization of the signal, so although Virgo signal-to-noise (SNR) is not as high as that of the LIGO detectors it's inclusion in the analysis still boosts sky localization. This was crucial for locating the position of the binary neutron star merger GW170817.
This method is effected by systematic uncertainties in the template waveform (i.e. the waveform that's compared to the observed signal) due to the breakdown of analytic approximations or numerical inaccuracies in simulating the waveform, and in calibration of the instruments which are independent at the different detectors, whereas the waveform errors are not. See section 5 of this.
In the Fisher Information Matrix Approximation (See section II.C), which can be accomplished with a single detector, the likelihood of the source parameters is approximated as a Gaussian centered on the true parameters, which is reasonable in the high SNR limit. This method is not so great for LIGO/Virgo because the SNR of detected events is not high (i.e., currently rarely above $\sim 20$). More generally, a Fisher matrix approach has a long history for source parameter estimation of coalescing binaries, see this seminal work by Cutler and Flanagan from 1994. The Fisher matrix is inverted to compute the covariance matrix which is used to compute the sky area. The full parameter space of a binary black hole is typically composed of 15 to 17 parameters, but only 9 parameters are needed for the approximate sky location.
For LISA specifically, Curt Cutler wrote a seminal paper about its expected angular resolution of a source's location in 1998 using the approximate Fisher matrix method. It's important to keep in mind that LISA works differently than LIGO, which is like a giant table-top interferometer where the lasers are set to precisely interfere destructively and when that interference is disturbed then a signal is detected. LISA will be in space, so it's not possible to precisely tune the lasers of the interferometer arms like is done with LIGO. Instead, elaborate and sophisticated methods have to be employed, see e.g. this. Anyway, as Cutler explains, "Information about the source position will be encoded in the measured
signal in three ways: (1) through the relative amplitudes and phases of the two polarization components, (2) through the periodic Doppler shift imposed on the signal by the detector’s motion around the Sun, and (3) through the further modulation of the signal caused by the detector’s time-varying orientation." For example, for supermassive black-hole binary mergers, he finds that LISA should achieve an angular resolution of (very roughly) $\sim$ 0.3 square degrees. Considering the many approximations he made in his analysis, preliminary work has shown that Cutler's predictions are remarkably accurate in most of the parameter space compared to a full Bayesian Fisher analysis. However, this is still probably not sufficient to precisely locate the host galaxy or cluster of the supermassive binary, since one square degree contains (at least) $10^4$ galaxies. Cutler argues that due to the large SNR of LISA, this might be sufficient however if many electromagnetic telescopes search within this square degree and might get lucky to see a flare from the merger if there happens to be enough gas in the vicinity of the merger.
In the future, a network composed of the LISA and Taiji detectors could utilize the timing triangulation method.