Video reconstruction of variable VLBI observations with neural fields

Supermassive black holes at the centre of active galactic nuclei drive the ejection of jets^1,2, highly relativistic, collimated streams of plasma that propagate through the intergalactic medium¹¹. Very long baseline interferometry (VLBI) is a primary tool to observe jets, enabling the detection of broad, unresolved features, called components, moving downstream at apparent superluminal speeds^12,13,14,15. Superluminal components have been associated with travelling shock waves that compress the magnetized plasma¹⁶, whereas alternative models explain the components with plasma instabilities and geometrical effects^17,18.

Traditional VLBI imaging relies on deconvolution techniques such as the widely adopted CLEAN⁵ algorithm. In recent years, alternative forward-modelling imaging approaches have been proposed^6,7,8, which achieve improved angular resolution by incorporating prior assumptions about the image in a parametrized model of the sky brightness¹⁸.

We present a forward-modelling VLBI imaging method, kine, that recovers a full-polarimetric video of variable radio sources. kine reconstructs the video by jointly recovering all frames simultaneously, learning and enforcing the spatio-temporal correlations in the data through a neural representation of the video, which infers the extent of the correlations during the imaging process. This method outputs a continuous representation of the brightness density distribution, sampleable at any time coordinate. By imaging multiple datasets simultaneously, each frame is reconstructed using substantially more information than in traditional, single-epoch imaging, leading to improved angular resolution and dynamic range that can overcome limited (u, v)-coverage or lower signal-to-noise ratio (SNR) on individual observations.

The continuity and high-resolution of our video allow us to apply post-processing techniques, such as optical flow^19,20, to recover the local instantaneous velocity field associated with the plasma. This is an advance over traditional kinematic analyses based on Gaussian model fitting^21,22 or wavelet component decomposition²³ because it resolves in space and time the projected velocity field describing the apparent motion of the emission of the source, rather than inferring motion by tracking discrete components. Assuming that spatial changes in the emission correspond to motion of the jet plasma, our method measures both the plasma bulk motion and the velocity of individual bright components, whereas previous methods could only measure the latter.

We showcase the performance of kine on 116 observations of blazar 3C 345 at 15 GHz, conducted from 1995 to 2022, by the Very Long Baseline Array (VLBA) as part of the MOJAVE monitoring program¹⁰. At a redshift of z = 0.593 and an inclination of 3°–6.8° from the line of sight^22,24, 3C 345 has been monitored for five decades^25,26, with recent observations showing a variable relativistic jet, characterized by the ejection of numerous superluminal components interpreted as moving shocks^{21,22,24,27,28,29}. With our imaging method, paired with motion analysis tools, we recovered a continuous, high-resolution, high-dynamic range polarimetric video of 3C 345 and measured the two-dimensional (2D) velocity field in the jet plasma, shedding light on the nature of jet dynamics and superluminal components.

The simultaneous imaging of all epochs with kine achieves a resolution of about 113 micro-arcseconds (μas), corresponding to a four-fold gain over the nominal value. The total achieved dynamic range is around 5 × 10⁵, about 140× more than traditional methods, which is due in part to the resolution improvement (approximately 14×) and in part to simultaneous imaging (about 10×). Figure 1 shows example frames out of the 116 video frames reconstructed (top rows). The full video is available in Extended Data Figs. 1 and 2 and Supplementary Information Video 1. The jet presents an inner jet, propagating westwards from a bright compact core, and a diffuse plume bending northwards, for a total extension of about 13 mas. Figure 1 (bottom rows) shows a magnification of the inner jet for three frames, with polarization ticks (left) and the plasma velocity vector field (right). From the jet core, plasma is continuously expelled downstream with occasional and noticeable ejection of discrete bright features, which initially move ballistically, before turning towards the jet axis, around 2–3 mas from the core, and vanishing into the plume (see, for example, the westernmost feature in Fig. 1b). These proper motions have been identified also in previous works through Gaussian model fitting at 15 GHz and other frequencies^21,29,30. We focus our analysis on the jet plasma kinematics as an example of the possible studies enabled by our method. Supplementary Information includes further analysis of the jet direction variability.

Jet velocity field

In contrast to the standard VLBI kinematic analysis, which relies on fitting Gaussian components to broad features travelling along the jet, we track the instantaneous, local velocity field of the plasma of the jet through optical flow estimation. This is enabled by the higher angular resolution of the reconstructions of our method, because of which we can obtain a velocity vector field, defined over space and time, at the same resolution as the video. Moreover, the time continuity of the video allows for motion-preserving interpolation between observed frames, circumventing the issue of irregular observation cadence. The instantaneous apparent velocity in the jet is shown in Fig. 1 (left), whereas in Fig. 2a, we present the time-average speed (i) and standard deviation (ii). By integrating the optical flow field from a given initial position, we obtain the trajectories followed by volume elements of the jet plasma, as shown in Fig. 2b. The optical flow trajectories following discrete components in the kine video agree with those obtained from model fitting of lower-resolution Gaussian components²¹ (Fig. 2c), which confirms the reliability of the optical flow method. However, the optical flow is not limited to tracking discrete bright components such as model fitting or wavelet decomposition. Owing to the high resolution of the video of kine and the locality of the optical flow vector field, it can also measure the velocity of any point in the jet, interpreted here as the velocity of the bulk jet flow in the absence of a passing component. Distinguishing between the speed of components and the speed of the flow is not possible with either Gaussian fitting or wavelet decomposition analysis. Details about frame interpolation with kine, a comparison between optical flow and model fitting and the optical flow limitations are discussed in the Supplementary Information.

From the optical flow measurements (Fig. 2d), the highest value of the mean apparent flow speed decreases from β_app = 12 ± 0.2c, between 1 mas and 3 mas from the core in the southern part of the jet, in which most of the components are ejected, to [9–11] ± 0.2c within the first 5 mas and [5–8] ± 0.2c further out in the diffuse emission. The maximum average speed corresponds to a physical speed of 0.997c for reported viewing angles. The speed standard deviation has no substantial location dependence and ranges between 3c and 6c, indicating that the plasma flow is turbulent and subject to relevant changes in apparent speed.

The velocities of the bright components (10–13c) are comparable with the average bulk plasma speed (9–12c) in the same region. This suggests that the components are not strongly shocked regions in the plasma flow, as previously proposed^16,31, but more likely localized over-densities in the plasma or regions with increased emissivity. In the presence of strong shocks, we would expect a significant difference between the bulk velocity of the shocked material and the pattern velocity of the shock, which is not the case for this source. Nevertheless, it is worth noting that similar values for the shock pattern speed and the plasma flow speed can also arise if protons are not highly relativistic³². We further test our hypothesis by analysing the polarimetric properties of the bright features.

Polarization and magnetic field

kine reconstructions yield videos of the Stokes parameters \({\mathcal{I}}\), \({\mathcal{Q}}\) and \({\mathcal{U}}\), showing the evolution of the linear polarization \(P={\mathcal{Q}}+i{\mathcal{U}}\) and fractional linear polarization \({m}_{{\ell }}:= | P| /{\mathcal{I}}\). Figure 1 shows the polarization field for three example frames, whereas Fig. 3a shows the time-averaged polarization field. An electric vector position angle (EVPA) structure in which the EVPAs align with the jet along the spine and orient perpendicularly near the jet border, recurs across multiple epochs (Fig. 1a) and is persistent also in the time-average polarization. When bright features are ejected, this pattern is disrupted as the feature travels through (Fig. 1b,c).

The EVPA structure is similar to the lower resolution polarization map reported in ref. ³³ and is consistent with a long-lasting toroidal magnetic field threading the jet in 3C 345 (refs. ^34,35,36), potentially resulting from a large-scale helical field dominated by its toroidal component³⁷.

Fractional polarization maps offer another probe into the nature of the travelling bright features. Figure 3b–e shows four such features, overlaid with instantaneous m_ℓ contours. If the features were propagating shocks, we would expect locally higher m_ℓ at their location¹⁶. Instead, although the passage of a bright feature produces a slight local disruption of the underlying polarization pattern, the peaks in m_ℓ show no spatial correlation with the features in total intensity.

The absence of this correlation further disfavours strong shocks in the jet. A travelling shock compresses the plasma, amplifying the magnetic field component parallel to the shock front because of the increased density. In a turbulent or toroidal magnetic field, the compression enhances field ordering and increases m_ℓ within the shocked region, whereas a predominantly poloidal field would produce the opposite effect^16,38. As our polarization pattern points to a toroidal magnetic field, if the bright components were shocks, we would expect localized increases in m_ℓ, yet this is not observed. Combining this with the pattern and bulk speed comparison, we interpret the bright, compact features as regions of enhanced emissivity associated with locally amplified magnetic fields, in which plasma turbulence increases the magnetic pressure and Doppler boosting further enhances their brightness when ejected southwards.

The analysis conducted in this work was enabled by two properties of kine. First, simultaneous imaging of all observations allows information to be shared across frames, improving resolution and dynamic range beyond what is achievable by frame-by-frame imaging (Supplementary Fig. 8). Second, the neural representation produces a smooth, continuous model of the flux density, sampleable at any time coordinate, enabling continuous and local motion analysis by optical flow. Applied to multi-epoch observations, kine can provide marked advances in jet kinematics studies: from high-resolution time-continuous videos, it becomes possible to measure the instantaneous local velocity field rather than only tracking model-fitted components. Future application to other blazar sources may lead to a reevaluation of the nature of assumed travelling shocks or to precise measurements of their pattern and flow speeds.

Multi-epoch blazar imaging is not the only possible application of kine. The algorithm was originally developed for horizon-scale single observations of Sagittarius A* with the Event Horizon Telescope (EHT). In that case, the source variability timescale is shorter than the full observation duration and the instantaneous coverage is too sparse to constrain individual frames, making dynamic imaging the only approach to resolve the variability of the source.

In static mode, kine achieves super-resolution comparable to recent regularized maximum likelihood and Bayesian methods (for example, eht-imaging⁷, Comrade⁸, Resolve⁶), improving over CLEAN⁵. In dynamic mode, unlike algorithms such as StarWarps³⁹, eht-imaging⁴⁰, ngMEM⁴¹, kine enforces correlations through implicit regularization rather than explicit morphological priors, making reconstructions robust and the algorithm source-agnostic. It also requires no careful tuning or expensive exploration of the regularization hyperparameters. Compared with dynamic Bayesian methods such as Resolve⁴², kine offers considerably shorter runtimes and better scalability to large datasets.

Generally, kine can image any VLBI dataset, single or multi-epoch, static or dynamic, in full polarization. It requires minimal expertise from the user, while implicit regularization limits human bias in the previous assumptions. Its ability to scale to large datasets with a short runtime makes it the only algorithm that can now reconstruct high-dynamic-range videos of multi-epoch observations. Beyond its potential for imaging dynamic sources, such as blazars, black holes or X-ray binaries, kine holds promise for VLBI tasks that require automation. Applications outside radio astronomy are also conceivable: MRI reconstruction poses formally the same problem as VLBI imaging, and kine could prove useful for imaging moving organs. Finally, the kine framework is naturally extensible to multi-frequency imaging by treating frequency as an additional input coordinate, making spectral index mapping a target for future development.

We presented a new dynamic imaging algorithm for VLBI observations of variable sources (Fig. 4) with an example application to multi-epoch observations of blazar 3C 345 from the MOJAVE program. The algorithm provides a time-continuous full-polarimetric video of the source at a resolution of about 113 μas and a dynamic range of approximately 5.1 × 10⁵. The super-resolution of kine is slightly higher than that of other algorithms, whereas the achieved dynamic range outperforms any available method because of the correlation of information across multiple observations. Combining our imaging results with optical flow results in a high-resolution map of the instantaneous projected velocity field in the source. For 3C 345, this allowed us to distinguish and compare the speed of discrete components and the average bulk flow speed.

We found no evidence that travelling bright components are strongly shocked regions, as previously proposed^16,31, because the component speeds are of the same order as the average flow speed in the same regions. The shock interpretation is also disfavoured by the absence of correlation between the bright features and the peaks in the fractional polarization. Therefore, we interpret the travelling features as Doppler-boosted regions of enhanced emissivity associated with local increases in magnetic pressure. This could not be concluded from earlier studies^21,29 because model fitting can only measure the speed of components, that is, the pattern speed of a possible shock, but not the underlying flow speed. Our results challenge only the current shock models for the 3C 345 blazar, but future systematic application of the proposed imaging and analysis technique to archival and future observations may produce more general results.

The kine algorithm is a static or dynamic image reconstruction algorithm for interferometric data. It is a forward imaging method, based on a neural field representation⁴³ of the brightness distribution of the source. In this section, we describe our methodology, starting with a general overview of VLBI measurements followed by an explanation of the kine imaging algorithm, which is shown in Fig. 4. Then, we describe the synthetic data used to test the performance of the algorithm and show the corresponding reconstructions. In the Supplementary Information, we present extensive tests validating the dynamic range, resolution, time interpolation and motion tracking of our method on synthetic datasets, of which we report the results here. We also present a performance comparison with the CLEAN method on real observations and the details of the optical flow velocity analysis.

VLBI measurements

In radio interferometric observations, each antenna in the array records a signal proportional to the received electromagnetic flux density. By the van Cittert–Zernike theorem, the time-averaged correlation product of the signals recorded by any pair of antennas (called visibility) is the Fourier transform of the flux density spatial distribution on the sky plane, evaluated at a frequency proportional to the distance of the antennas. For a pair of antennas A and B, with projected baseline vector \(\bar{b}=({b}_{x},{b}_{y})\), observing at wavelength λ, the ideal complex visibility V_AB(u, v, t) is related to the flux density distribution \({\mathcal{I}}(x,y,t)\) by

where \((u,v)=\left(\frac{{b}_{x}}{\lambda },\frac{{b}_{y}}{\lambda }\right)\) are the x and y spatial frequencies (refer to ref. ⁴⁴ for a detailed description of radio interferometry observations and imaging). The above equation, written for Stokes \({\mathcal{I}}\), holds for all Stokes visibilities \(({V}_{\mathrm{AB}}^{({\mathcal{I}})},{V}_{\mathrm{AB}}^{({\mathcal{Q}})},{V}_{\mathrm{AB}}^{({\mathcal{U}})},{V}_{\mathrm{AB}}^{({\mathcal{V}})})\).

In practice, different sources of noise corrupt the measurement of visibilities. They are classified as baseline-dependent errors and site-dependent errors, so the measured visibilities can be expressed as

where \({G}_{{\rm{A}},{\rm{B}}}{{\rm{e}}}^{{\rm{i}}{\phi }_{{\rm{A}},{\rm{B}}}}\) are the site-dependent errors in amplitude and phase, referred to as complex gains, and ϵ_AB is the thermal noise, which is Gaussian-distributed with baseline-dependent standard deviation⁴⁴. The most problematic source of error is the complex gains, as they might be incorrectly estimated from the a priori calibration, whereas thermal noise can be fully characterized and incorporated in the loss function through the uncertainties σ_AB of the visibilities.

Therefore, although complex visibilities are the fundamental data product resulting from VLBI observations, the imaging process can use different data products that are constructed to be unaffected by site-dependent amplitude or phase corruptions. kine supports the following data products: complex visibilities V_AB, amplitudes of the visibilities |V_AB|, closure phases \({\varPhi }_{\mathrm{ABC}}:= \arg ({V}_{\mathrm{AB}}{V}_{\mathrm{BC}}{V}_{\mathrm{CA}})\), closure amplitudes A_ABCD ≔ |V_ABV_CD|/|V_ACV_BD| and the logarithm of closure amplitudes.

Model

kine is a forward imaging algorithm for interferometric data, which models the polarized flux density distribution of a continuous image or video using a neural field⁹, parameterized by weights \({\mathcal{W}}\). The neural field is implemented as a coordinate-based neural network, specifically a multi-layer perceptron (MLP)⁴⁵, which takes spatial coordinates (right ascension x and declination y) and time (t) as input, and outputs the total intensity flux density \(\hat{I}\), fractional linear polarization \({\widehat{m}}_{{\ell }}\), electric vector position angle \(\widehat{\chi }\) and fractional circular polarization \({\hat{m}}_{c}\) at each sky location and time, as shown in Fig. 4:

The hat symbol denotes estimated quantities, as opposed to the true or observed properties of the source. The Stokes parameters are then computed from the polarization quantities \((\hat{{\mathcal{I}}},{\hat{m}}_{{\ell }},\hat{\chi },{\hat{m}}_{c})\) using the transformation:

We compute VLBI measurements from the predicted video using a fully differentiable forward model and optimize the MLP weights by minimizing a χ² data-fit loss between the predicted and observed measurements. Within this forward model, estimated complex visibilities \(\widehat{V}\) and associated data products are obtained by evaluating the 2D Fourier transform of the network-estimated flux density \(\hat{{\mathcal{I}}}\) at the (u, v, t) points corresponding to the observation tracks. In practice, \(\hat{{\mathcal{I}}}\) is evaluated on a discrete set of (x, y, t) points, and the visibilities are therefore computed using a discrete Fourier transform F implemented as a matrix multiplication. For each observation time t_j and u–v point i, the estimated visibilities are

where k runs over all spatial points at which the video is estimated, and both the estimated total intensity and the estimated visibilities are dependent on the network parameters \({\mathcal{W}}\).

MLP networks are known to exhibit a ‘spectral bias’⁴⁶, by which training with gradient descent preferentially captures low spatial frequencies. In the context of VLBI imaging, this behaviour is advantageous, as it acts as an implicit regularizer that promotes smooth reconstructions and suppresses spurious high-frequency structure in both space and time. Moreover, the continuous nature of the neural representation enables the video reconstruction to be evaluated at times when no observations are available, by leveraging correlations learnt from neighbouring frames. This implicit regularization through neural representations has been successfully applied to other ill-posed inverse problems in astrophysical imaging, in which many possible solutions are consistent with a given set of measurements^47,48,49.

Optimization

We formulate the optimization as a noise-weighted data-fitting problem. Assuming Gaussian noise residuals in the data products, which is valid for visibilities and amplitudes and for closure quantities in the high-SNR regime, the natural loss function \({\mathcal{L}}\) is given by the sum of the χ² contributions from the selected data products:

where D indicates a generic data product, the index j runs over all the observed times, from 1 to the total number of observed times N_t; the index i runs over all data corresponding to the observed time t_j, from 1 to the total number of data N_D,j; and k_D is a normalization factor that takes into account the degrees of freedom of the data product. Any combination of one or more kinds of data products can be used in the loss function, and the loss function terms associated with each data product are weighted equally.

The optimization process consists of iterating the following steps until convergence. First, the algorithm estimates the video, with the current network parameters \({{\mathcal{W}}}_{k}\), on a coordinate grid with x and y regularly spaced and t spaced according to the observation times. Then it computes the loss function from the estimated video and the observed data. Finally, it updates the MLP parameters to \({{\mathcal{W}}}_{k+1}\), by applying a gradient descent optimization step of the loss function. More specifically, we use the Adamax optimizer⁵⁰. In this work, convergence is considered to be reached when the moving average of the loss function over a window of 100 iterations does not decrease by more than 1% for at least 1,000 iterations, provided that all data products in the loss function have reached a value of about 1. This convergence was achieved after approximately 2 × 10⁴ iterations, a value that we finally set as the total number of iterations.

Calibration-free imaging details

When imaging a dataset with no gain corruption, the best data products to use are complex visibilities, because they contain the full amount of available information. In that case, the imaging procedure in kine consists of running the optimization algorithm estimating all Stokes parameters simultaneously. However, if the dataset contains amplitude and/or phase gains, imaging must be performed using closure quantities. When imaging with closure or log-closure amplitudes, the information about the total flux of the frames is lost so a total flux regularizer term is added to the loss function, constraining the total flux in each frame to a value either provided by the user or computed from the instantaneous visibility of the shortest baseline. When imaging with closure phases, the information about the absolute location of the source in the frame is lost, so the reconstructed movie may show the source drifting smoothly across the frame. We solve this issue post-imaging by re-aligning the frames to a chosen feature, such as a jet core. The alignment procedure is detailed in the Supplementary Information.

The current version of kine does not model station gains simultaneously with the imaging. For the total intensity reconstruction, this problem is solved by relying on closure quantities. The closure quantities for the \({\mathcal{Q}}\) and \({\mathcal{U}}\) signals, however, are gain independent only in the absence of polarization leakage⁵¹, but, even in that case, the low SNR of the \({\mathcal{Q}}\) and \({\mathcal{U}}\) data, combined with the loss of information due to the use of closure quantities, may result in insufficient data constraints. Therefore, we found that the most effective procedure to produce a full-polarimetric video, in the presence of gain corruption, is to first image the data in total intensity only, using the appropriate closure quantities, then self-calibrate the data to the total intensity video reconstruction, and finally, image the self-calibrated data in full Stokes, using complex visibilities.

The kine imaging method was developed for video reconstruction of a variable source, but it can reduce to static imaging by removing the time dimension from the input coordinates.

Architecture

kine is implemented in Python using the JAX deep learning framework⁵². We use a four-layer MLP for total intensity imaging and a deeper six-layer MLP for full-Stokes imaging of the self-calibrated data, both with 256 nodes per layer and a residual skip connection⁵³ between the first and the last layers. The node outputs in all layers except the last are normalized with BatchNorm⁵⁴ across the full set of input coordinates to have zero mean and unit standard deviation. Hidden layers use the Gaussian Error Linear Unit activation⁵⁵, whereas the output layer uses SoftPlus for \(\hat{{\mathcal{I}}}\) and a sigmoid for \({\widehat{m}}_{{\ell }}\), \({\widehat{m}}_{c}\) and \(\widehat{\chi }\). Standard interferometry functions are handled by the eht-imaging library⁷. During optimization, spatial input coordinates are sampled on a regular grid, whereas the time coordinates follow the observation timestamps. For the imaging of 3C 345, we used 200 × 200 spatial input points with 75 μas spacing.

The runtime depends on the number of visibilities and epochs, the field-of-view to beam-size ratio, and the data products used. Optimizing the network on 116 epochs of 3C 345 data from the VLBA took approximately 1.3 h on four NVIDIA A100 GPUs, whereas imaging 4 days of M 87* EHT observations took around 20 min on a single A100.

Initialization

When the (u, v)-coverage of the VLBI data is sufficiently dense, such as for VLBA observations, kine is able to reconstruct the image or video starting from a random initialization of the network weights. However, if imaging with closure phases, it may be useful to initialize the network to a simple shape, such as a Gaussian, to constrain the image to the correct location in the frame. We found that, although it is important that the initialization image has the correct total flux and covers approximately the area in which most of the flux of the source should be located, the specific shape of the initialization image does not affect the resulting image or video obtained when convergence is reached.

The initialization is carried out by optimizing the network directly on the initialization image \({{\mathcal{I}}}_{\mathrm{init}}(x,y)\), with a pixel-to-pixel distance loss function:

without the computation of the image Fourier transform or other data products. After the weights of the network are optimized to output the initialization image, the optimization on the data products begins.

Synthetic data validation

We validated the algorithm on multiple synthetic datasets prepared to resemble the morphology, coverage and noise of the real 3C 345 data. The first ground truth video consists of a geometric jet model in which the emission direction of the jet is precessing along a circular orbit, emitting ballistic components, which expand and decrease in brightness as they move away from the core. This results in a model with a continuous winding structure, which we refer to as model1. The second model, model2A, is constructed in a similar way, but with longer spacing between the emitted components, so that the model is composed of discrete, separate components emitted at different angles. The components have an initial full width at half maximum (FWHM) = 80 μas, which increases slightly as they move away from the core at a constant apparent speed of β = 11.4. The third model, model3A, is similar to the second but with the addition of a parabolic jet profile as a background to simulate a more realistic jet structure. We used these three models for the main validation tests, but we also prepared variations of model2A and model3A, with components of different sizes—FWHM = 60 μas for model2B and model3B and FWHM = 160 μas for model2C—and moving at different speeds—β = 22.7 for model2D and β = 5.7 for model2E—to test the resolution and motion recovery. For each observed epoch, we simulated a synthetic dataset from the model, using the corresponding real (u, v)-coverage and applying the same amount of thermal noise. We did not need to introduce complex gain corruption, because the closure quantities used in the first step of the imaging process are invariant under gain errors.

We imaged the synthetic data with the same parameters used for real data. Here, as an example, we present and discuss the reconstructions of model1, model2A and model3A. The ground truth videos are shown in the first, third and fifth rows of Extended Data Fig. 3, whereas the reconstructions of the synthetic data are in the second, fourth and sixth rows of the same figure. These show that kine is able to recover the models with high fidelity, which is reflected by the high values of the average normalized cross-correlation between the ground truth frames and the reconstructed ones (Extended Data Table 1).

In the Supplementary Information, we present the details of additional validation tests that quantify the dynamic range, super-resolution, component characterization, time interpolation and motion tracking achieved by our method, and we report the quantitative results of the tests in the subsequent sections. The improvements in reconstruction quality achieved by kine over traditional imaging methods, especially the dynamic range and resolution improvements, are dependent on the quality and quantity of the observations, so they should not be intended as intrinsic performance gains independent of dataset quality and temporal coverage.

Forward-modelling imaging algorithms can achieve higher resolution than traditional inverse methods by incorporating prior information. From the validation tests, we see that kine dynamic reconstructions have effective resolutions between about 100 μas and 125 μas, with an average of 113 μas, whereas kine static reconstructions have effective resolutions between about 110 μas and 145 μas, with an average of 125 μas. This shows that the main resolution improvement is achieved already at the static reconstruction level, with an improvement factor of 3.8 over the nominal resolution of 475 μas, whereas the dynamic reconstruction provides a marginally higher improvement factor of 4.2.

The dynamic range improvement of kine over traditional methods is a result of the combination of two factors: the resolution improvement of kine and the simultaneous imaging of multiple epochs. From the validation tests, we conclude that kine static reconstructions achieve a dynamic range between 6 × 10³ and 8 × 10⁴, whereas dynamic ones achieve a dynamic range from 2 × 10⁵ to 7 × 10⁵. This shows that dynamic imaging provides a higher dynamic range than static imaging. In this case, the improvement factor is of the order of about 10, which is compatible with the \(\sqrt{N}=\sqrt{116}\simeq 11\) limit.

We tested the ability of kine to correctly resolve and locate small discrete components in the core region by fitting a double 2D Gaussian model on kine reconstructions. The tests show that components can be precisely located starting at about 110 μas from the core and tracked throughout their whole path. For a component fully detached from the core, the position angle can be recovered within 1° of the true value, whereas for those connected to the core, it is recovered correctly within 7°.

We tested the limit of the velocity that could be recovered with kine combined with the optical flow, and we conclude that we can accurately measure apparent velocities in the 3C 345 dataset up to at least approximately 23c, which exceeds the maximum speeds observed in the real data.

Owing to the continuous nature of neural field representations, kine reconstructions can be sampled at any arbitrary time, independent of the observed epochs. For data with the coverage and quality of 3C 345 observations, kine provides good, motion-preserving time interpolation for frames up to 6 months apart from the nearest observations. In the present work, we use time-interpolated frames in the optical flow analysis to have a video sampled at regular intervals. The longest separation between an interpolated frame and the nearest observation is 5.7 months, which is within the range of reliable interpolation.

Comparison of CLEAN and kine

As shown from the validation tests, even when the instantaneous coverage is sufficient for static imaging, dynamic reconstructions achieve higher resolution and dynamic range than static ones. Here, we estimate the effective resolution and dynamic range of the static and dynamic reconstructions of kine of the real 3C 345 observations and compare them with images obtained with CLEAN. Supplementary Fig. 8 shows an example frame comparing traditional static imaging with CLEAN, snapshot static imaging with kine and dynamic imaging with kine. Both kine reconstructions show consistency with the CLEAN image, although they are able to resolve considerably more structure in the jet. As we do not have a ground truth in this case, we can estimate the dynamic range only with the traditional definition, while we estimate the resolution using the power spectrum comparison with the CLEAN model convolved with the nominal beam.

We obtain an average dynamic range of 3.6 × 10³ for CLEAN, 4.9 × 10⁴ for kine static and 5.1 × 10⁵ for kine dynamic, resulting in an improvement of one order of magnitude for kine static over the traditional CLEAN imaging method and two orders of magnitude for kine dynamic. The dynamic ranges of kine static have a wide distribution of values, whereas that of kine dynamic varies less. This means that kine dynamic is effective in improving the reconstruction quality of epochs with lower data quality by leveraging the information from previous and subsequent epochs.

We obtain a resolution of 141 μas for kine static and a resolution of 106 μas for kine dynamic, representing an improvement factor of 3.4 for kine static and 4.5 for kine dynamic with respect to CLEAN, whose average nominal resolution of 475 μas is set by the diffraction limit. These results are consistent with the estimates from the tests on synthetic data, although we consider the resolution estimate from real data slightly less reliable than that from synthetic data because the comparison is not carried out with the ground truth but rather with the CLEAN reconstruction. Therefore, we prefer to rely on the more conservative estimate from synthetic data for our conclusions.

The observed data used in the imaging of 3C 345 are publicly available on the MOJAVE website (https://www.cv.nrao.edu/MOJAVE/index.html). The kine video reconstruction, the synthetic data, the post-imaging analysis and plotting notebooks are available at GitHub (https://github.com/mariannafoschi/kine3c345/).

The specific version of the kine algorithm used in this work is publicly available at GitHub (https://github.com/mariannafoschi/kine3c345/).

1. Blandford, R. D. & Znajek, R. L. Electromagnetic extraction of energy from Kerr black holes. Mon. Not. R. Astron. Soc. 179, 433–456 https://doi.org/10.1093/mnras/179.3.433 (1977).

   Article  ADS  Google Scholar 

2. Blandford, R. D. & Payne, D. G. Hydromagnetic flows from accretion discs and the production of radio jets. Mon. Not. R. Astron. Soc. 199, 883–903 https://doi.org/10.1093/mnras/199.4.883 (1982).

   Article  ADS  Google Scholar 

3. Narayan, R., Igumenshchev, I. V. & Abramowicz, M. A. Magnetically arrested disk: an energetically efficient accretion flow. Publ. Astron. Soc. Jpn. 55, L69–L72 (2003).

   Article  ADS  Google Scholar 

4. Tchekhovskoy, A., Narayan, R. & McKinney, J. C. Efficient generation of jets from magnetically arrested accretion on a rapidly spinning black hole. Mon. Not. R. Astron. Soc. Lett. 418, L79–L83 (2011).

   Article  ADS  Google Scholar 

5. Högbom, J. A. Aperture synthesis with a non-regular distribution of interferometer baselines. Astron. Astrophys. Suppl. Ser. 15, 417–426 (1974).

   ADS  Google Scholar 

6. Junklewitz, H., Bell, M. R., Selig, M. & Enßlin, T. A. Resolve: a new algorithm for aperture synthesis imaging of extended emission in radio astronomy. Astron. Astrophys. 586, A76 (2016).

   Article  ADS  Google Scholar 

7. Chael, A. A. et al. Interferometric imaging directly with closure phases and closure amplitudes. Astrophys. J. 857, 23 (2018).

   Article  ADS  Google Scholar 

8. Tiede, P. Comrade: composable modeling of radio emission. J. Open Source Softw. 7, 4457 (2022).

   Article  ADS  Google Scholar 

9. Xie, Y. et al. Neural fields in visual computing and beyond. Comput. Graph. Forum 41, 641–676 (2022).

   Article  Google Scholar 

10. Lister, M. L. et al. MOJAVE. XV. VLBA 15 GHz total intensity and polarization maps of 437 parsec-scale agn jets from 1996 to 2017. Astrophys. J. Suppl. Ser. 234, 12 (2018).

    Article  ADS  Google Scholar 

11. Oei, M. S. S. L. et al. Black hole jets on the scale of the cosmic web. Nature 633, 537–541 (2024).

    Article  ADS  CAS  PubMed  Google Scholar 

12. Cohen, M. H. et al. Radio sources with superluminal velocities. Nature 268, 405–409 (1977).

    Article  ADS  Google Scholar 

13. Vermeulen, R. C. & Cohen, M. H. Superluminal motion statistics and cosmology. Astrophys. J. 430, 467–494 (1994).

    Article  ADS  Google Scholar 

14. Gómez, J.-L., Marscher, A. P., Alberdi, A., Jorstad, S. G. & García-Miró, C. Flashing superluminal components in the jet of the radio galaxy 3C120. Science 289, 2317–2320 (2000).

    Article  ADS  PubMed  Google Scholar 

15. Jorstad, S. G. et al. Polarimetric observations of 15 active galactic nuclei at high frequencies: jet kinematics from bimonthly monitoring with the very long baseline array. Astron. J. 130, 1418–1465 (2005).

    Article  ADS  Google Scholar 

16. Marscher, A. P. & Gear, W. K. Models for high-frequency radio outbursts in extragalactic sources, with application to the early 1983 millimeter-to-infrared flare of 3C 273. Astrophys. J. 298, 114–127 (1985).

    Article  ADS  CAS  Google Scholar 

17. Raiteri, C. M. et al. Blazar spectral variability as explained by a twisted inhomogeneous jet. Nature 552, 374–377 (2017).

    Article  ADS  CAS  PubMed  Google Scholar 

18. Fuentes, A. et al. Filamentary structures as the origin of blazar jet radio variability. Nat. Astron. 7, 1359–1367 (2023).

    Article  ADS  Google Scholar 

19. Horn, B. K. P. & Schunck, B. G. Determining optical flow. Artif. Intel. 17, 185–203 (1981).

    Article  Google Scholar 

20. Meinhardt-Llopis, E., Pérez, S. J. & Kondermann, D. Horn–Schunck optical flow with a multi-scale strategy. Image Processing On Line 3, 151–172 (2013).

    Article  Google Scholar 

21. Lister, M. L. et al. Monitoring of jets in active galactic nuclei with VLBA experiments. XVIII. Kinematics and inner jet evolution of bright radio-loud active galaxies. Astrophys. J. 923, 30 (2021).

    Article  ADS  CAS  Google Scholar 

22. Weaver, Z. R. et al. Kinematics of parsec-scale jets of gamma-ray blazars at 43 GHz during 10 yr of the VLBA-BU-BLAZAR program. Astrophys. J. Suppl. Ser. 846, 98 (2022).

    Google Scholar 

23. Mertens, F. & Lobanov, A. Wavelet-based decomposition and analysis of structural patterns in astronomical images. Astron. Astrophys. 574, A67 (2015).

    Article  ADS  Google Scholar 

24. Steffen, W., Zensus, J. A., Krichbaum, T. P., Witzel, A. & Qian, S. J. A helical model for the compact jet in 3C345. Astron. Astrophys. 302, 335–342 (1995).

    ADS  Google Scholar 

25. Unwin, S. C. et al. Superluminal motion in the quasar 3C 345. Astrophys. J. 271, 536–550 (1983).

    Article  ADS  Google Scholar 

26. Biretta, J. A., Moore, R. L. & Cohen, M. H. The evolution of the compact radio source in 3C 345. I. VLBI observations. Astrophys. J. 308, 93–109 (1986).

    Article  ADS  CAS  Google Scholar 

27. Lister, M. L. et al. MOJAVE. XVII. Jet kinematics and parent population properties of relativistically beamed radio-loud blazars. Astrophys. J. 874, 43 (2019).

    Article  ADS  CAS  Google Scholar 

28. Pötzl, F. M. et al. Probing the innermost regions of AGN jets and their magnetic fields with RadioAstron. IV. The quasar 3C 345 at 18 cm: magnetic field structure and brightness temperature. Astron. Astrophys. 648, A82 (2021).

    Article  Google Scholar 

29. Röder, J., Ros, E., Schinzel, F. K. & Lobanov, A. P. Up around the bend: a multiwavelength view of the quasar 3C 345. Astron. Astrophys. 684, A211 (2024).

    Article  ADS  Google Scholar 

30. Jorstad, S. G. et al. Kinematics of parsec-scale jets of gamma-ray blazars at 43 GHz within the VLBA-BU-BLAZAR program. Astrophys. J. 846, 98 (2017).

    Article  ADS  Google Scholar 

31. Blandford, R. D. & Königl, A. Relativistic jets as compact radio sources. Astrophys. J. 232, 34–48 (1979).

    Article  ADS  CAS  Google Scholar 

32. Marscher, A. P. Relativistic jets in active galactic nuclei. AIP Conf. Proc. 856, 1–22 (2006).

    Article  ADS  Google Scholar 

33. Pushkarev, A. B. et al. MOJAVE - XX. Persistent linear polarization structure in parsec-scale AGN jets. Mon. Not. R. Astron. Soc. 520, 6053–6069 (2023).

    Article  ADS  CAS  Google Scholar 

34. Lyutikov, M., Pariev, V. I. & Gabuzda, D. C. Polarization and structure of relativistic parsec-scale AGN jets. Mon. Not. R. Astron. Soc. 360, 869–891 (2005).

    Article  ADS  Google Scholar 

35. Fuentes, A., Gómez, J. L., Martí, J. M. & Perucho, M. Total and linearly polarized synchrotron emission from overpressured magnetized relativistic jets. Astrophys. J. 860, 121 (2018).

    Article  ADS  Google Scholar 

36. Fuentes, A., Torregrosa, I., Martí, J. M., Gómez, J. L. & Perucho, M. Magnetized relativistic jets and helical magnetic fields. II. Radiation. Astron. Astrophys. 650, A61 (2021).

    Article  ADS  Google Scholar 

37. Gabuzda, D. C., Reichstein, A. R. & O’Neill, E. L. Are spine-sheath polarization structures in the jets of active galactic nuclei associated with helical magnetic fields? Mon. Not. R. Astron. Soc. 444, 172–184 (2014).

    Article  ADS  Google Scholar 

38. Hughes, P. A., Aller, H. D. & Aller, M. F. Polarized radio outbursts in BL Lacertae. II. The flux and polarization of a piston-driven shock. Astrophys. J. 298, 301–315 (1985).

    Article  ADS  CAS  Google Scholar 

39. Bouman, K. L. et al. Reconstructing video of time-varying sources from radio interferometric measurements. IEEE Trans. Comput. Imaging 4, 512–527 (2018).

    Article  Google Scholar 

40. Johnson, M. D. et al. Dynamical imaging with interferometry. Astrophys. J. 850, 172 (2017).

    Article  ADS  Google Scholar 

41. Mus, A. & Martí-Vidal, I. New-generation maximum entropy method: a Lagrangian-based algorithm for dynamic reconstruction of interferometric data. Mon. Not. R. Astron. Soc. 528, 5537–5557 (2024).

    Article  ADS  Google Scholar 

42. Arras, P. et al. Variable structures in M87* from space, time and frequency resolved interferometry. Nat. Astron. 6, 259–269 (2022).

    Article  ADS  Google Scholar 

43. Mildenhall, B. et al. NeRF: representing scenes as neural radiance fields for view synthesis. In Computer Vision – ECCV 2020 (eds Vedaldi, A., Bischof, H., Brox, T. & Frahm, J.-M.) 405–421 (2020).

44. Thompson, A. R., Moran, J. M. & Swenson, G. W. Jr. Interferometry and Synthesis in Radio Astronomy 3rd edn (Springer, 2017).

45. Popescu, M.-C., Balas, V. E., Perescu-Popescu, L. & Mastorakis, N. Multilayer perceptron and neural networks. WSEAS Trans. Circuits Syst. 8, 579–588 (2009).

    Google Scholar 

46. Rahaman, N. et al. On the spectral bias of neural networks. In Proc. 36th International Conference on Machine Learning Vol. 97, 5301–5310 (PMLR, 2019).

47. Levis, A., Chael, A. A., Bouman, K. L., Wielgus, M. & Srinivasan, P. P. Orbital polarimetric tomography of a flare near the sagittarius A* supermassive black hole. Nat. Astron. 8, 765–773 (2024).

    Article  ADS  Google Scholar 

48. Zhao, B., Levis, A., Connor, L., Srinivasan, P. P. & Bouman, K. L. Single view refractive index tomography with neural fields. In Proc. 2024 IEEE/CVF Conference on Computer Vision and Pattern Recognition 25358–25367 (IEEE, 2024).

49. Zhong, E. D., Bepler, T., Berger, B. & Davis, J. H. CryoDRGN: reconstruction of heterogeneous cryo-em structures using neural networks. Nat. Methods 18, 176–185 (2021).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

50. Kingma, D. & Ba, J. Adam: A method for stochastic optimization. In Proc. International Conference on Learning Representations (ICLR, 2015).

51. Roelofs, F. et al. Polarimetric geometric modeling for mm-vlbi observations of black holes. Astrophys. J. Lett. 957, L21 (2023).

    Article  ADS  Google Scholar 

52. Bradbury, J., Frostig, R., Hawkins, P., Johnson, M. J. & Katariya, Y. JAX: composable transformations of Python+NumPy programs (v.0.3.13). GitHub. http://github.com/jax-ml/jax (2018).

53. He, K., Zhang, X., Ren, S. & Sun, J. Deep residual learning for image recognition. In Proc. 2016 IEEE Conference on Computer Vision and Pattern Recognition 770-778 (IEE, 2016).

54. Ioffe, S. & Szegedy, C. Batch normalization: accelerating deep network training by reducing internal covariate shift. In Proc. 32nd International Conference on International Conference on Machine Learning—Volume 37 (eds Bach, F. & Blei, D.) 448–456 (JMLR, 2015).

55. Hendrycks, D. & Gimpel, K. Gaussian Error Linear Units (GELUs). Preprint at arxiv.org/abs/1606.08415 (2016).

Download references

Author notes

1. These authors contributed equally: Marianna Foschi, Brandon Zhao, Antonio Fuentes

Authors and Affiliations

1. Instituto de Astrofísica de Andalucía (IAA-CSIC), Granada, Spain

   Marianna Foschi, Antonio Fuentes & José L. Gómez

2. Division of Physics, Mathematics and Astronomy, Caltech, Pasadena, CA, USA

   Marianna Foschi

3. Department of Computing and Mathematical Sciences, Caltech, Pasadena, CA, USA

   Brandon Zhao & Katherine L. Bouman

4. Department of Computer Science, University of Toronto, Toronto, Ontario, Canada

   Aviad Levis

5. David A. Dunlap Department of Astronomy & Astrophysics, University of Toronto, Toronto, Ontario, Canada

   Aviad Levis

Contributions

M.F., B.Z. and A.F. developed and coded the imaging method, applied it to the observations, performed all quantitative analysis and wrote the manuscript. K.L.B. and A.L. conceived the method and provided guidance on the algorithm and code development and results discussion. J.L.G. guided the overall structure and scientific content of the manuscript, supervised the analysis and physical interpretation of the results, and provided guidance on the algorithm development.

Corresponding author

Correspondence to Marianna Foschi.

Competing interests

The authors declare no competing interests.

Peer review information

Nature thanks Yosuke Mizuno and the other, anonymous, reviewer(s) for their contribution to the peer review of this work.

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

Reprints and permissions