Time-dependent visibility modelling of a relativistic jet in the X-ray binary MAXI J1803$-$298

In an interferometer, each pair of telescopes measures a complex visibility, which, according to the van Cittert-Zernike theorem (van Cittert, 1934; Zernike, 1938) is an element of the Fourier transform of the sky brightness distribution. Typical VLBI observations can be up to several hours in length, so that as the Earth rotates, the orientations of the separation vectors (the baselines) of the pairs of telescopes projected onto the plane of the sky change. This samples more unique visibilities, increasing the overall sensitivity and quality of the image reconstruction.

Since the complex visibility plane (often called the $uv$-plane) can never be completely sampled, the inverse Fourier transform of the visibilities is a convolution of the true sky brightness distribution and the inverse Fourier transform of the sampling function. Many imaging algorithms have been developed with the aim of reconstructing the true sky brightness distribution from the incomplete information. In radio astronomy, the standard technique is the CLEAN algorithm (Högbom, 1974; Schwarz, 1978; Clark, 1980), which is a deconvolution technique that represents the sky as a sum of point sources, and attempts to iteratively subtract out the artefacts and side lobes of the sampling function from the inverse Fourier transform of the visibilities.

There is another class of techniques that first attempt to reconstruct a version of the sky brightness distribution, before comparing that reconstruction to the underlying data. Examples of this class of methods are the so called maximum entropy methods (MEM) or regularised maximum likelihood (RML) methods, that try to solve for the best possible image by fitting the image pixels to the data while providing constraints via the use of regularization terms that favour certain features in the image, e.g. entropy, sparsity, or smoothness (Frieden, 1972; Cornwell & Evans, 1985; Narayan & Nityananda, 1986). These methods are not as popular as the CLEAN algorithm, although in recent years they have gained attention, particular by groups such as the Event Horizon Telescope (EHT) collaboration (e.g. Chael et al., 2016, 2018; Akiyama et al., 2017; EHT Collaboration et al., 2019; Broderick et al., 2020; Broderick et al., 2022).

Another example of this approach is model fitting. In this technique, simple model source components with analytic representations in the Fourier domain (e.g. point sources or Gaussians) are fit directly to the visibilities (e.g. Shepherd et al., 1994; Martí-Vidal et al., 2014), greatly reducing the number of free parameters in the imaging problem. Historically, model fitting was first used before imaging, with early two element interferometers (e.g. Fomalont, 1968). Model fitting has been used to study transient jets launched by LMXBs, since they are often seen in images as compact point sources or Gaussians (e.g. Miller-Jones et al., 2019).

One of the fundamental assumptions of VLBI is that over the length of an observation the target source is non-variable. Jets launched by LMXBs can travel across the resolution element of an image in a matter of minutes (e.g. Wood et al., 2021) and can vary by a significant fraction of their flux density on the same time-scale (e.g. Miller-Jones et al., 2019), violating this assumption. The simplest solution is time binning (e.g. Fomalont et al., 2001; Miller-Jones et al., 2019), where the full observation is split into short time bins, within which the source is relatively static, each to be imaged individually. This technique requires the source to be bright enough so that it can be significantly detected in each time bin, since within a single time bin the sensitivity and $uv$-coverage are greatly reduced, making this technique difficult for standard LMXB ejecta, which are typically only tens of mJy in brightness.

Recently, more sophisticated techniques have been developed that seek to improve upon the time binning procedure. In Wood et al. (2021) we described a dynamic phase centre tracking technique by which we applied an incremental phase shift to each time bin of an observation before stacking the time bins back together to effectively ’track’ a jet component with a given proper motion. Other recent developments have been focused on capturing the variability of the super-massive black holes M87^∗ and Sagittarius A^∗ in EHT observations (e.g. Bouman et al., 2017; Johnson et al., 2017; Arras et al., 2022). These approaches aim to extend MEM and RML methods to simultaneously reconstruct images from all of the time bins in an observation, while explicitly enforcing continuity across the full set of images, to leverage information from the entire observation to enhance the quality and sensitivity of each individual image. The further development of techniques that can capture intra-observational variability is key to making more precise measurements of the proper motions and ejection dates of transient jets, and thus determining the causal connection between changes in the inner accretion flow and the launching of relativistic jets in LMXBs.

MAXI J1803$-$298 (hereafter J1803) was first discovered as a new X-ray transient in the early stages of an outburst on 2021 May 1 (MJD 59335) by the Monitor of All-sky X-ray Image (MAXI; Serino et al., 2021) nova alert system. It was quickly localised by both NICER (Gendreau et al., 2021a, b) and Swift, the latter of which also detected an optical counterpart (Gropp et al., 2021). Spectroscopy with the Southern African Large Telescope suggested that J1803 was an LMXB (Buckley et al., 2021). NICER, AstroSat, and NuSTAR X-ray spectral analysis further suggested that J1803 was an accreting stellar-mass black hole, as opposed to an accreting neutron star (Bult et al., 2021; Jana et al., 2021; Xu & Harrison, 2021), viewed relatively edge-on, with an inclination above $70\degree$. On 2021 May 4 (MJD 59338.9), Espinasse et al. (2021) first detected J1803 at radio wavelengths with MeerKAT. On 2021 May 11 (MJD 59345), 10 days after its initial discovery, AstroSat detected a state transition of J1803 (Jana et al., 2021), with their observations suggesting that the source had entered the hard-intermediate state. MAXI/GSC was unable to observe J1803 for eight days from 2021 May 4 (MJD 59338), but on 2021 May 12 (MJD 59346) they also reported that the source was in the intermediate state, with the transition to the soft state occurring on 2021 May 28 (MJD 59362; Shidatsu et al., 2022). J1803 remained in the soft state for $\sim 5$ months, with the reverse soft-to-hard state transition occurring between 2021 October 13 and 2021 October 19 (MJD 59500-59506; Steiner et al., 2021).

We present the results of a radio monitoring campaign of J1803 during the state transition, with the Very Long Baseline Array (VLBA), the Australia Telescope Compact Array (ATCA), and the Atacama Large Millimeter/Sub-Millimeter Array (ALMA). In order to account for the intra-observational variability of our VLBA observations, we have developed a new model-fitting approach in which we jointly fit a single time-evolving model to all of the visibilities in a single observation, rather than on a time bin by time bin basis. This allowed us to leverage all of the information from a full observation in a single fit to constrain the motion and flux density variability of the detected components. We first demonstrate our validation of this technique with synthetic data sets designed to replicate the typical variability we would expect in our observations, before presenting the results of our application of this technique to our VLBA observations of J1803.

The paper is organised as follows. We describe our observations, calibration, imaging, and model fitting procedure in Section 2. We present the results of this analysis in Section 3. We discuss our results in Section 4 and present our conclusions in Section 5.

We first imaged our VLBA data within aips using the CLEAN algorithm. On the days when there were multiple observations, we imaged each of these epochs separately, before concatenating them to increase sensitivity and $uv$-coverage, since they were separated by a short amount of time. We refer to these concatenated epochs by their first letter (e.g. epochs B1 and B2 became epoch B). J1803 is close ($\sim 4\degree$) to the Galactic centre, and thus the longer baselines were affected by scattering due to the dense, turbulent interstellar medium, resulting in angular broadening. In order to recover images with a resolution that matched the effective resolution of the angularly broadened data, we applied a Gaussian $uv$-taper with 30% power at 50 mega-wavelengths (the typical maximum baseline of these observations was 150–250 mega-wavelengths). While we tried many uv-taper sizes, we chose this uv-taper to maximize the recovered flux density in the images by suppressing the scattered long baselines, while not compromising the image quality by removing too many inner baselines. Thus, the resolution, as marked in the lower left corner of the final images (Fig. 1), is lower than in typical VLBA observations. We observed J1803 using the best fitting position from the Swift/XRT localisation (Gropp et al., 2021), although we shifted the phase centre of all of the observations to align with the centroid position of the compact component detected in epoch A. We were only able to reliably detect J1803 in epochs A, B, and C, likely the result of the transient jet ejecta adiabatically expanding and fading as they became resolved out. We henceforth focus on these three observations. We were unable to perform any self calibration to improve the significance of our detections, since the source was too faint in all of our VLBA epochs.

Following our basic characterisation of the images in aips, we successfully fit elliptical Gaussian model components directly to the complex visibilities in DIFMAP for each of the three VLBA epochs. The signal to noise was too low in each of these epochs to split the individual observations into time bins to track the intra-observational variability of the detected components. While model fitting is able to reduce the number of free parameters in the imaging problem relative to algorithms like CLEAN, this approach is still limited by having to perform model fitting on multiple individual time bins in order to detect motion and flux density variability. To overcome this issue, and to constrain the nature of the intra-observational variability in our VLBA data, we have developed a new model fitting approach in which we fit a single time-evolving model to the full set of visibilities from an observation.

By parameterising the motion and flux density variability of modelled emission components, we are able to leverage information from all of the time bins in a single fit, rather than having to individually fit distinct models to each of the separate time bins. Since transient jets often appear in images as compact point sources or Gaussians, we are able to use these simple emission profile models, which have analytical representations in the visibility domain. We therefore fit these simple models directly to the time-stamped interferometric visibilities. We can allow any of the parameters that describe these model components (e.g. position, size, flux density) to be time-dependent, and explicitly parameterise their variability with analytic expressions that are included in the model. We can then predict the time-varying visibilities of our model and compare these to the measured visibilities of our observations, to fit these time-variable models directly to the underlying data without requiring any Fourier transforms.

In this work, we fit the position of the source components (modelled as circular Gaussians) with a ballistic (i.e constant velocity) model, in a polar coordinate system. While later-time deceleration of jet ejecta has been seen (e.g. Espinasse et al., 2020), on these scales, a constant velocity model is adequate. Our equations of motion for the position of the source at time $t$ were therefore,

and

where $\Delta x(t)$ and $\Delta y(t)$ are the positions of the component, relative to the phase centre of the observations, in the directions of Right Ascension (RA) and Declination (Dec.), respectively; $x_{0}$ and $y_{0}$ are the fit positions of the source at the beginning of the observation, relative to the phase centre; $t_{0}$ is the time at the beginning of the observation; $\dot{r}$ is the proper motion of the fitted source, $\theta$ is the position angle on the sky along which the model component moves in °East of North; and $\delta$ is the declination of the source. The phase centre of the observations was shifted to align with the centroid position of the component detected in the image of epoch A. In our model, we allowed the flux density to vary linearly as,

where $F_{0}$ is the flux density of the source at the beginning of each observation and $\dot{F}$ is its derivative. We can easily make this model more sophisticated by adding further degrees of freedom to account for variability in other model parameters, for example by allowing for expansion. However, in this work we only consider simple models that evolve linearly with time in position and flux density, due to the low signal to noise and particularly sparse uv-coverage of our observations. Therefore our fit parameters for this model (and associated units) are $x_{0}$ (mas), $y_{0}$ (mas), $\dot{r}$ (mas hr^-1), $\theta$ (deg.), $F_{0}$ (mJy), $\dot{F}$ (mJy hr^-1), and the full-width half maximum angular size of the circular Gaussian component, $\theta_{\text{FWHM}}$ (mas).

In order to perform parameter estimation for these time-evolving models, we used the framework of Bayesian inference (see van de Schoot et al., 2021, for a primer on Bayesian inference). Due to the phase referencing, we assumed that the data were well calibrated and thus were not corrupted by incorrect station-based complex gains, and instead only considered thermal noise. Since the thermal errors on complex visibilities are Gaussian (Thompson et al., 2017), we used a Gaussian likelihood in our parameter estimation. We describe our model formalism and application of Bayesian inference in the visibility domain in Appendix A. Station-based complex gains could be included in the model as fitted parameters, or could even be avoided completely by the use of calibration-independent closure quantities as the data products (e.g. Chael et al., 2018); however we do not consider that in this work. In order to explore the posterior probability distribution to estimate the best fitting model parameters, we used the Bayesian inference algorithm nested sampling (Skilling, 2006) implemented in the dynesty³³3https://github.com/joshspeagle/dynesty (Speagle, 2020) Python package. Nested sampling is well suited for our model fitting requirements, given its ability to efficiently traverse multi-modal posterior distributions. We used the python library eht-imaging⁴⁴4https://achael.github.io/eht-imaging/ (Chael et al., 2018) for the handling of our UVFITS data and for pre-processing. Our implementation is available via GitHub⁵⁵5https://github.com/Callan612/MAXIJ1803-Model-Fitting.

We also analyze available X-ray data on J1803 to track the evolution of the outburst, as an auxiliary source of information guiding the interpretation our VLBA observations. This includes publicly available light curves from MAXI/GSC⁶⁶6http://maxi.riken.jp/ (Matsuoka et al., 2009), and Swift/BAT⁷⁷7https://swift.gsfc.nasa.gov/results/transients/ (Krimm et al., 2013), along with data from Swift/XRT.

We reduced Swift/XRT data and extracted a light curve in the 0.3-10 keV band for J1803 using HEASOFT (v6.29) and the UK Swift Science Data Centre online platform⁸⁸8https://www.swift.ac.uk/user_objects/ (Evans et al., 2007, 2009). For extraction, we used the source coordinates from Gropp et al. (2021) and binned the light curve by averaging per each individual observation. Hardness was calculated as the ratio between count rate in the 1.5-10 keV band over that in the 0.3-10 keV band.

We describe two of our synthetic observations to demonstrate the capability of our model fitting technique. We used model parameters that are similar to the parameters that we found best described our real observations of J1803. First we demonstrate the model fitting results for a circular Gaussian moving with a proper motion of 1 mas hr^-1, along the $y$-axis ($\theta=0\degree$). Fig. 5 shows the joint and marginal posterior probability distributions for the fit to this synthetic data set with the true model parameters overlaid. In all but two parameters, the model fitting is able to recover the true parameter values within a $1\sigma$ credible interval, and in all of them within $2\sigma$.

We also demonstrate the model fitting results for an identical circular Gaussian, but with no motion. Fig. 6 shows the joint and marginal posterior probability distributions for this fit. The key identifier for a lack of motion in the source is in the marginal posterior for $\dot{r}$, which is consistent with having a mode at the boundary value of zero. In the case where there is no motion, the position angle $\theta$ along which the component moves should be uniformly distributed. In the high signal-to-noise case we found this to be true, however there is a clear structure in the marginal posterior distribution for $\theta$ as seen in Fig. 6. This is likely the result of the sparse $uv$-coverage and poor signal to noise of the simulated observation. Motion on the plane of the sky manifests in the data as a change in the slope of the phase of the complex visibilities. When the signal is faint and the $uv$-coverage is sparse, as is the case with this synthetic observation, it is harder to constrain motion along the direction in which the $uv$-coverage is most sparse, resulting in the structure in the posterior distribution for $\theta$ seen in Fig. 6. This is analogous to imaging, where the shape and orientation of the resolution element of the image (which is described by the restoring beam) is determined by the distribution of the $uv$-coverage, which results in poorer resolution in the direction in which the $uv$-coverage is most sparse. Synthetic data sets generated with more complete $uv$-coverage and higher signal to noise showed a marginal probability distribution for $\theta$ that was much more uniform.

We present, in Fig. 1, our CLEAN images and a visualisation of our model fitting results for epochs A, B, and C. In Table 2 we detail the results of our model fits. We discuss the identification and behaviour of these VLBA components in Sections 4.2 and 4.3. The image of Epoch A consisted of a single compact source with a peak intensity of $1.2\pm 0.1$ mJy/beam, where the uncertainty is the $1\sigma$ statistical uncertainty reported by the AIPS task JMFIT. In epoch A we allowed for a moving circular Gaussian, however we found that the marginal posterior probability distribution for $\dot{r}$ was consistent with having a mode at zero, with a $1\sigma$ credible upper limit on the motion of the component in epoch A being $\dot{r}<0.48$ mas hr^-1. The marginal posterior probability distribution for $\theta$ (Fig. 7) is not uniform and shows two peaks, the narrowest and largest at $\sim 150\degree$ east of north and a small but broader peak at $\sim 10\degree$ east of north. The posterior probability in between these peaks is non-negligible. The bi-modal peaks of this posterior distribution approximately correspond to the position angle of the CLEAN synthesised beam, i.e. the directions along which the $uv$-coverage is more sparse. This is identical to the behaviour seen in the synthetic observation described in Section 2.3.1, suggesting that this component is most likely stationary. The best fit location of this component was,

While we found the component in epoch A was most likely stationary, it was rapidly rising in flux density, with $F_{0}=0.76_{-0.18}^{+0.19}$ mJy and $\dot{F}=0.59\pm 0.12$ mJy hr^-1. The circular Gaussian component had a full-width half-maximum of $5.2_{-0.4}^{+0.5}$ mas. We tried splitting epoch A in half, and performed model fitting on each half independently, and found that we were able to consistently constrain the lack of motion and the rapid rise in flux density of the source in both halves of the observation. This epoch was the faintest of the three epochs, however it enjoyed the most sensitivity due to its observation duration.

In the image of epoch B, we found a single component $\sim 10$ mas to the south-east (to the left in the rotated image) of the component in epoch A. It was also slightly extended with an asymmetric flux density distribution skewing towards the south-east. The peak intensity of the component was $4.0\pm 0.2$ mJy/beam. In epoch B we found evidence of motion in our model fitting, with $\dot{r}=1.37\pm 0.14$ mas hr^-1 at a position angle of $168\pm 4\degree$ east of north. The position angle along which the components in epochs A, B, and C lie is approximately $135\degree$ east of north. The full-width half-maximum of the Gaussian component in epoch B was $4.4\pm 0.2$ mas. Epoch B had the longest lever arm in time to detect motion, since in epoch A we observed geodetic blocks at the beginning and end of the observation. Epoch B, however, enjoyed less sensitivity than epoch A since it was comprised of two $\sim 1$ hour long epochs separated by a $\sim 1.5$ hour gap.

The image of Epoch C consisted of a single component $118$ mas away from the component in epoch A, to the north-west (right in the rotated image). Similar to epoch B, the component also appeared to be slightly resolved, although not as extended as the source in epoch B. It had a peak intensity of $3.8\pm 0.3$ mJy/beam. Since we were unable to adequately converge on a solution for a moving circular Gaussian component, we fit a single static circular Gaussian with a full-width half-maximum of $5.5\pm 0.2$ mas, at a separation of $122.0_{-1.7}^{+1.9}$ mas from the component in epoch A at a position angle of $313.9_{-0.6}^{+0.8}\degree$ east of north. As marked in the images, the best model fit position of the source is $\sim 4$ mas further away from component A than the component in the image. It is unclear why this is the case, although it could be because of some faint, extended, asymmetric emission which we cannot detect in the imaging but can constrain in the model fitting. The image of epoch C contains some bright fringes with a direction and angular separation consistent with originating from the shortest baseline (LA-PT), possibly hinting at the existence of some larger extended emission that is only detected on the shortest baseline. The image of epoch D was dominated by this fringing, and we were unable to detect any compact source structure, which may be the result of the component in epoch C expanding and only being detected on this shortest baseline. In epoch C, the circular Gaussian component had a full-width half-maximum size of $5.5\pm 0.2$ mas. The model fits for epoch B and epoch C were both consistent with having flat light curves, unlike epoch A where we saw a rapid rise in the flux density of the source over the length of the observation. Epoch C was the shortest and least sensitive of the three epochs in which J1803 was clearly detected, which could explain why we were unable to constrain a time-evolving model with motion.

We produced images for epochs A, B, C, and D, having flagged out the LA-PT baseline, which removed the large scale fringing. We found that for epoch C, the integrated flux density in the image was slightly reduced. We were still unable to detect any compact components in epoch D. Our modelling results were completely consistent with, and without the LA-PT baseline.

By stacking all of the observations of J1803 in the soft state (MJD 59362–59500; Shidatsu et al., 2022), we were able to put a 3$\sigma$ upper limit on the flux density of the core in the soft state of $<0.097$ mJy (at 8.2 GHz). In the final observation, epoch P, we were unable to detect the core of the system in the hard state with a 3$\sigma$ upper limit of $0.105$ mJy (at 4.9 GHz).

In Fig. 2 we show the radio flare at the peak of the outburst, observed in both our ATCA and VLBA observations. Since we fit a linearly evolving flux density model to our VLBA observations, we plot the flux density at the beginning and end of each of the three VLBA epochs in which the source was detected. In the same figure we also show the evolution of the spectral index, $\alpha$ ($S_{\nu}\propto\nu^{\alpha})$, over the ATCA observations. In the first ATCA observation at the beginning of the radio flare, the source spectrum was flat ($\alpha=-0.1\pm 0.1$). In the second ATCA observation, J1803 reached its peak flux density at both 5.5 and 9 GHz, and had the steepest spectrum, with $\alpha=-0.8\pm 0.2$. The source then began to fade, with the spectral index flattening gradually but remaining negative. In the first ATCA observation, there was no evidence of intra-observational variability over the length of the observation. In the second observation, which lasted for $\sim 1.5$ hours, both the 5.5 and 9 GHz flux densities were rising steadily, from 8.6/5.5 mJy to 10.0/6.3 mJy (5.5/9 GHz), suggesting that this observation did not correspond to the true peak of the radio flare. We also include the ALMA epoch 2, 98.5  GHz flux density measurement, which occurred very close to the VLBA epoch C observation. In Fig. 4, we show the image of the ALMA epoch 2, with the positions of the components from the VLBA epochs A, B, and C marked on top. The position of the unresolved component detected by the ALMA is consistent with the position of the component detected in the VLBA epoch C. We split the ALMA epoch 2 into 30 second time-bins, which revealed the source had an approximately constant flux density over the length of this observation, which was consistent with the flat flux density profile of the nearby VLBA epoch C.

The VLBA observations provided denser sampling of the flux density of J1803 around the peak of this flare than the ATCA observations, with epoch A showing rapid intra-observation brightening at 8.4 GHz, just prior to the peak of the ATCA radio flare. Our model fitting revealed no significant intra-observational flux density evolution in VLBA epochs B and C. J1803 was in the intermediate state for the entirety of the time-span shown in Fig. 2.

In Fig. 3 we show the full ATCA light curves with the ALMA and VLBA measurements, accompanied by the Swift/BAT, Swift/XRT, and MAXI/GSC 1-day averaged light curves of J1803 around the peak of the outburst and the state transition. In the ATCA light curves we observed a re-brightening following the initial radio flare, beginning at $\sim$ MJD $59365$ and peaking on MJD $59384.8$. During this period of re-brightening, we did not detect any emission from a compact core or transient radio jet in our VLBA observations.

Our model fitting approach extends the traditional model fitting implementations of software like DIFMAP (Shepherd et al., 1994) or UVMULTIFIT (Martí-Vidal et al., 2014) by parameterising the variability and motion of model components over a full observation. In this way we were able to leverage information from a full observation to constrain this variability. Model fitting approaches have captured the variability of jets launched by LMXBs, for example in V404 Cygni. However, Miller-Jones et al. (2019) performed their model fitting separately on each individual time bin. This was possible for V404 Cygni because all of the individual components were bright ($10^{1}$-$10^{3}$ mJy) point-like sources. Our new approach is similar in premise to the new time-resolved imaging techniques developed to create movies from EHT observations of the super-massive black holes M87^∗ and Sgr A^∗, which seek to leverage information from the full observation to enhance the quality of each time-binned image by enforcing or parameterising continuity between time bins (Bouman et al., 2017; Johnson et al., 2017; Arras et al., 2022). Rather than trying to perform a full pixel-by-pixel image reconstruction, as these techniques do, we parameterise our data with simple model components, which greatly reduces the number of free parameters in the reconstruction. The use of these simple models is physically motivated, since transient jets launched by LMXBs often appear as point sources or compact Gaussian components in VLBI images. We therefore do not require the more sophisticated image reconstructions that have been developed to account for the complex turbulent flows and asymmetry in objects such as M87^∗ and Sgr A^∗.

We validated our technique using synthetic observations that were designed to replicate the typical source behaviour and variability we see in observations of jets from LMXBs. We found that we were able to reliably recover the input model parameters for a range of different synthetic data sets.

With our initial CLEAN imaging we were able to measure the locations, sizes and flux densities of discrete components in our observations. However we were unable to characterise the variability of these components, since the signal-to-noise of the observations was too low, and the $uv$-coverage was too sparse for time binning, even with standard model fitting within DIFMAP. With our new model fitting approach, we have been able to measure the proper motions of the components in epochs A and B by fitting a ballistic proper motion model, as well as measuring the flux density evolution of the components in epochs A, B, and C, by fitting a linear model. We were unable to fit more complicated models, including models with acceleration/deceleration, non-linear flux density evolution, elliptical Gaussian components, or expanding Gaussian components, due to the sparse $uv$-coverage and low signal-to-noise of the observations. In observations with more complete $uv$-coverage and higher signal-to-noise, we could in future extend our technique to include these more sophisticated models. For consecutive observations where the same components are detected multiple times, we could even use this technique to perform a single fit across multiple observations. In the case of J1803, the measured positions and proper motions (or lack thereof) for the different components in the three VLBA observations in which we have a detection suggest that they are distinct.

From our VLBI imaging and model fitting, we have identified three separate components in epochs A, B, and C, which we will refer to as components 1, 2, and 3, respectively. Based on our model fitting and imaging we identify the three components as follows.

Our modelling suggested that component 1 is most likely stationary, and is rapidly rising in flux density. Given its apparent lack of motion, its compact structure, its presence towards the beginning of the outburst, and the fact that component 2 is approximately moving away from its position, we suggest that this component is likely the core of the system. It is not clear if the radio emission is originating from a compact, steady jet prior to it becoming quenched in subsequent epochs, or if the rapidly rising flux density originates from the rise of radio emission from slow-moving transient ejecta. While we suspect that component 1 is the core, we were unable to confirm the position of the core with any other observations. In epoch P, we were unable to detect the core of the system in the hard state, nor were we able to detect the core by stacking together all of the epochs in which there were non-detections. No other instrument that observed J1803 provides enough angular resolution to give an independent confirmation of the location of the core. Our only constraint on the position of the core of the system is therefore our measured position in epoch A, given in Section 3.

Component 2 is a jet moving away from the location of the core, which we believe to be component 1. We note that the direction of the motion of component 2 does not point directly back to component 1. It is not clear if this is due to some misalignment between the images, some directional bias in the modelling introduced by the non-uniform nature of the $uv$-coverage, or if component 1 is truly not the core. We also note that the markers plotted on top of the images in Fig. 1 do not capture the uncertainty of the model fits. If component 1 is not the core, then our only constraint on its position is that it should be along the axis of the motion of component 2, in the opposite direction.

Although we were unable to constrain the motion of component 3 (or lack thereof), we suggest it is most likely a jet, given its transient nature, and its position $\sim 120$ mas away from component 1 approximately opposite from the direction in which component 2 moves. Our failure to constrain models that parameterised the motion of component 3 is likely the result of the low signal-to-noise ratio of the observation as well as the fact that the jet was starting to become resolved out during this epoch.

Since we do not detect the core in epoch B, and assuming that component 1 is the core of J1803, we can constrain the time period during which the core switched off to be between MJD 59347.49 and MJD 59348.34. We constrained the ejection date of component 2 to be within this period of time, however it is not clear if the quenching of the compact jet in the core coincided with the ejection of component 2. In H1743-322, Miller-Jones et al. (2012) observed that the quenching of the compact radio jet emission occurred during or immediately after the ejection of transient ejecta. In J1803, following the ejection of component 2, the compact radio core was quenched for the remainder of the intermediate state and into the soft state. We were unable to detect the reestablishment of the radio core in our VLBA observations.

In epoch A, we observed the core to be rapidly rising in flux density at the beginning of the radio flare prior to the ejection of component 2. As discussed in Section 4.3.1, the rise of these radio flares is often attributed to either the expansion of a jet or the acceleration of jet material as it moves away from the launch site or the formation of internal shocks, and thus the rise generally is expected to lag the ejection of the jet material. The rapid rise in epoch A could be due to a sudden brightening of the compact jet prior to the ejection of component 2, the adiabatic expansion of a much slower moving optically thick jet component close to the core, or internal shocks in jet material deposited close to the core of J1803. If component 1 is not the core, its rapid rise could be attributed to the evolution of an isolated, slow moving jet component, which may be seen later as either component 2 or 3. We note that this component would have to undergo some acceleration prior to epoch B for this to correspond to component 2. We cannot rule out that this may be a slow-moving component 3 present in epoch A, since we have no constraints on the motion of component 3.