© The Authors 2024

1. Introduction

On the 9 October 2022, all satellites equipped for transient detection were triggered by the extraordinary γ-ray burst (GRB) 221009A (Veres et al. 2022; Bissaldi et al. 2022; Ursi et al. 2022; Piano et al. 2022; Gotz et al. 2022; Frederiks et al. 2022; Tan et al. 2022; Mitchell et al. 2022; Liu et al. 2022; Lapshov et al. 2022; Xiao et al. 2022; Ripa et al. 2022; Dichiara et al. 2022; Kennea et al. 2022). At a redshift of z = 0.151 (de Ugarte Postigo et al. 2022; Malesani et al. 2024), GRB 221009A holds the record of the highest ever measured isotropic equivalent energy (E_γ, iso ≳ 10⁵⁵ erg – Lesage et al. 2023). It is the brightest GRB in the last 50 years and it is estimated to be a one in ∼10 000 years occurrence based on the observed flux distribution of other known long GRBs (O’Connor et al. 2023; Burns et al. 2023; Malesani et al. 2024). Such a unique event initiated an unprecedented follow-up campaign, characterised by extensive temporal and spectral coverage. At the highest energies, the LHAASO Collaboration reported the detection of sustained emission well above 1 TeV (LHAASO Collaboration 2023; Cao et al. 2023). At the lower end of the electromagnetic spectrum, radio observations of GRB 221009A commenced just three hours post-burst and detected the brightest ever radio counterpart, reaching a flux density of 60 mJy (Bright et al. 2023). Initial attempts to model the multi-wavelength afterglow emission considered contributions from both the reverse shock (RS) and the forward shock (FS) resulting from the deceleration of the ultra-relativistic jet by the surrounding material (Ren et al. 2023; Sato et al. 2023; Laskar et al. 2023; O’Connor et al. 2023; Gill & Granot 2023; Zheng et al. 2024). However, uncertainties persist in the final interpretation of the data, despite incorporating most of the presently known physical ingredients governing the dynamics and emission of GRB jets.

Unique measurements able to independently constrain the afterglow evolution can be obtained with milliarscsecond resolution observations. Very long baseline interferometry (VLBI) allows for direct measurements of the size of the emission region, together with high-precision astrometry. As a result, proper motion and source expansion can be measured (Taylor et al. 2004; Mooley et al. 2018; Ghirlanda et al. 2019). If the viewing angle θ_v between the observer line of sight and the GRB jet axis is smaller than the jet half-opening angle θ_j (‘on-axis’ GRB), the projected image during the afterglow is expected to expand, but not to show appreciable proper motion. Conversely, if the outflow is observed ‘off-axis’ (θ_v > θ_j), an apparent superluminal motion is expected. To date, measurements of the size and expansion of the emitting region have only been possible for GRB 030329 (Taylor et al. 2004, 2005), providing the first direct evidence of the relativistic expansion of GRB outflows. Over the last two decades, numerous campaigns were aimed to repeat the success of GRB 030329 (e.g. Nappo et al. 2017; Salafia et al. 2022; Giarratana et al. 2022). However, no event shone brightly and long enough to allow for an expansion measurement. On the other hand, for the multi-messenger event GW 170817 (Abbott et al. 2017a,b; Margutti & Chornock 2021), VLBI observations were fundamental to measure the apparent superluminal motion and to constrain the size of the emitting region of the non-thermal electromagnetic counterpart (Mooley et al. 2018; Ghirlanda et al. 2019), proving, for the first time ever, that the mergers of two neutron stars are able to successfully launch ultra-relativistic jets.

Here we present our VLBI follow-up campaign on GRB 221009A. The data reduction is detailed in Sec. 2. The method implemented to measure the source properties from radio observations is described in Sec. 3. In Sec. 4 we present the results of our campaign and discuss the physical implications in Sec. 5. Throughout the work, we assume Planck Collaboration VI (2020) cosmological parameters. With these parameters, the angular diameter distance at z = 0.151 is d_A = 560.3 Mpc. Therefore, 1 mas separation corresponds to 2.72 pc in projection at such a distance.

2. VLBI observations and data reduction

2.1. European VLBI Network

We observed the field of GRB 221009A with the European VLBI Network (EVN) from 40 to 261 days post-burst (PI: Giarratana, project code: RG013). Given the target-of-opportunity nature of the proposal, not all antennas were available at all epochs. Table A.1 lists the antennas joining each epoch. Table 1 presents a summary of the properties of the observations. The observations were performed in two different bands centred at 4.9 and 8.3 GHz. The data were recorded at 4 Gbits s⁻¹. Dual polarisation products (RR, LL) were correlated at the Joint Institute for VLBI in Europe (JIVE, Dwingeloo, Netherlands) using the Super FX Correlator (SFXC; Keimpema et al. 2015) into sixteen sub-bands with 32 MHz bandwidth and 64 channels each. For the last epoch, RG013 F, the data were correlated into eight sub-bands with 32 MHz bandwidth and 64 channels each. The first, EVN epoch (RG013 A), carried out 6 days post burst at a central frequency of 22.2 GHz, was not usable due to unfavourable observing conditions.

The observations consisted of phase-referencing cycles with 4.5 and 2.5 minutes on the target at 4.9 and 8.3 GHz, respectively, and 1.5 minutes on the phase calibrator. Further scans every approximately 30 minutes on some ‘check’ sources were also included. Throughout the observations, some scans on a fringe finder were performed. The radio source J190536.4+194308 (J1905+1943 hereafter) and the Very Large Array Sky Survey (VLASS) compact radio source J191142.50+195200 (J191142+1952 hereafter) were used as phase calibrators in the first two (RG013 B and C) and in the last three observations (RG013 D, E and F), respectively.

The calibration was performed using AIPS¹ (Greisen 2003), following the standard procedure for EVN phase-referenced observations². The amplitude calibration, which accounts for the bandpass response, the antenna gain curves and the system temperatures, was performed by applying the gains derived by the EVN pipeline. We performed a correction for the dispersive delay using the IONEX files from the International GNSS Service (vlbatecr procedure in AIPS), we calculated a manual single band delay on the fringe finder (vlbampcl procedure in AIPS) and we carried out the global fringe fitting on the phase calibrator (fring task in AIPS) using a model of the source derived by a concatenation (in CASA, McMullin et al. 2007) and self-calibration (in Difmap, Shepherd et al. 1994) of all the visibilities on the source obtained across the various epochs. Solutions were interpolated (clcal task in AIPS) and applied to the phase calibrator itself, the target and some check sources (see Appendix A). For the last three epochs, we corrected the visibilities of J191142+1952 by fixing the phase centre in CASA to the actual position of the phase calibrator, as the initial position of this phase calibrator was not constrained with a sub-mas resolution.

Images of the sources were produced using Difmap. For the analysis presented in this paper, we selected the image with the best signal-to-noise ratio (S/N) among the two images produced before and after the self-calibration of the phase calibrator, respectively. For RG013 C, the flux density of the GRB enabled a self-calibration in phase in AIPS (solint = 1 min). Further information on the structure and the data reduction process can be found in Appendix A.

2.2. Very Long Baseline Array

The Very Long Baseline Array (VLBA) data were acquired between 44 and 262 days post-burst (PI: Atri, project code: BA160). The central frequency was 15.2 GHz, with a total bandwidth of 512 MHz, divided into 4 spectral windows of 128 MHz and 256 channels each, in full circular polarisation (RR, LL, RL, LR). The number of participating stations contributing useful data was 7, 8, 10 and 10 in experiments BA160 B, C, C1 and D respectively (see Table A.1). Each observation included approximately 30-minute-long geodetic-style blocks at the beginning and at the end of the observation, used to determine troposphere modelling errors. The central part of the observations included scans on fringe finder bright calibrators and repetitions of a J1905+1943 – J1925+2106 – GRB 221009A sequence, where J1905+1943 and J1925+2106 are known VLBA calibrators, with respective durations of 30s – 30s – 80s.

The data were correlated at the NRAO in Socorro using the Distributed FX software correlator (DiFX; Deller et al. 2011). The data reduction was carried out in AIPS, following standard procedures for continuum phase-referencing experiments. Procedures vlbaeops, vlbaccor, vlbampcl, vlbabpss, vlbaamp were carried out in this order for the initial bandpass and amplitude calibration. The following step consisted in the calibration of the troposphere modeling errors by running the task fring on the geodetic blocks, followed by mbdly and delzn. The final phase, rate, and delay fringe-fitting was carried out separately on J1905+1943 and J1925+2106, yielding high S/N and well-behaved solutions for both sources. The solutions from the closer phase calibrator, J1905+1943, were applied to the target field. After preparing a model of the phase calibrator using Difmap, a cycle of amplitude and phase solutions were determined for the calibrator itself and applied to the target to further refine the calibration. Finally, we produced single-source frequency-averaged datasets for the target, which were imaged in AIPS with a natural weighting scheme.

Our VLBA campaign included one more epoch, BA160 A, at approximately 14 days post-burst. However, as the antennas were pointed at an incorrect position in the sky, the GRB fell outside the primary beam of the VLBA, which is approximately 3 arcmin at 15 GHz. While the reduced sensitivity (approximately 25% of nominal) still allowed for the detection of the burst, a satisfactory calibration of the complex visibilities was hampered. Therefore, we did not include this experiment in our analysis.

3. Methods

3.1. Source flux density, size and average apparent expansion velocity estimate

In order to extract information about the total flux density, size and position of the source from each of our epochs, we fitted a circular Gaussian source model to the calibrated visibility data adopting a Markov Chain Monte Carlo (MCMC) approach. The method, which is an extension of that adopted in Salafia et al. (2022), is described in Appendix B. The projected angular diameter of the source image is proportional to the full width at half maximum (FWHM) of the fitted circular Gaussian, with a proportionality constant of order unity that depends on the detailed surface brightness profile (Granot et al. 1999, 2005; Taylor et al. 2004; Pihlström et al. 2007; Granot 2008; Salafia et al. 2022). In what follows, we set the proportionality constant equal to 1, and discuss it whenever relevant. Once this diameter is measured, the average apparent expansion velocity can be calculated (assuming the size to be zero at the time t₀ of the explosion) as

$$\begin{array}{c}\hfill ⟨{\beta }_{\mathrm{app}}⟩=\frac{(1+z)d_{\mathrm{A}}s}{2(t_{\mathrm{obs}}-t_0)c},\end{array}$$(1)

where s is the FWHM, t_obs is the time of the observation, and c is the speed of light.

Table 1 summarises the result of the circular Gaussian fitting, along with the derived average apparent expansion velocity. In Appendix B we provide more detailed information in the form of corner plots that illustrate the posterior probability density of the flux density and source size from the circular Gaussian fitting. Figure 1 additionally shows ‘violin plots’ that illustrate a kernel density estimate of the posterior probability density of the FWHM for each epoch.

Fig. 1. Source size as a function of time. The source size constraints obtained as described in Section 3.1 are shown in the form of violin plots of different colours, centred at the observing time of each epoch and of proportional width to the posterior probability density of the FWHM. In addition, we show the median and 68% credible interval with an error bar of the same color, or the 95% credible upper limit with a triangle if the former interval extends to 0. The black dashed line and the two grey shaded areas show respectively the median, 68% credible interval and 95% credible interval of the posterior predictive distribution of the source size evolution obtained from fitting a power law model s ∝ tobsa to the sizes from all the epochs. The inset shows the posterior probability density of the slope a from such fit.

3.2. Source size evolution model fitting

In order to fit a size evolution model s_m(t_obs, θ) to the observations, where θ is a vector of free parameters, we adopted a Bayesian approach. By Bayes’ theorem, and given the fact that the size estimates from different observations are independent, the posterior probability on θ is proportional to the prior π(θ) times the product of the likelihoods. This can be written as

$$\begin{array}{c}\hfill P\left(\mathit{\theta }|\{{\mathit{d}}_i{\}}_{i=1}^M\right)\propto \pi (\mathit{\theta }){\displaystyle \prod _{i=1}^M}\frac{P\left(s_{\mathrm{m}}(t_{\mathrm{obs},i})|{\mathit{d}}_i\right)}{\pi (s)},\end{array}$$(2)

where M is the number of epochs included in the fit, d_i is the data (i.e. the visibilities) of the i-th epoch, t_obs, i is the time of the i-th observation, π(s) = Θ(s) (where Θ is the Heaviside step function) is the prior on the size adopted in the circular Gaussian fits, P(s | d_i) is the posterior from such fits (Eq. B.3) marginalised on all parameters except s. In order to evaluate Eq. (2), we approximated the marginalised posterior on the size P(s | d_i) with a Gaussian kernel density estimate based on the posterior samples derived from the MCMC described in Sect. 3.1. This allowed us to sample the posterior on θ again through an MCMC approach.

3.3. Forward and reverse shock size evolution and proper motion model

In order to interpret our observations in the context of the standard afterglow scenario, we derived a simple physical model of the size evolution and, in the case of a jet not aligned with the observer’s line of sight, the proper motion of the source expected if the emission is dominated by either the FS or the RS produced as a relativistic jet expands into an external medium with a power law number density profile n(R) = A(R/R_⋆)^−k, where R is the distance from the explosion site (i.e. the progenitor vestige) and R_⋆ = 5.5 × 10¹⁷ cm is a reference radius³. We assumed a uniform jet angular energy profile for simplicity, with an isotropic-equivalent kinetic energy E, a half-opening angle θ_j, an initial Lorentz factor Γ₀ and a duration T (which sets the jet radial width ΔR ∼ cT). The viewing angle is assumed to be either θ_v = 0 (on-axis, for the calculation of the projected size) or θ_v > θ_j (off-axis, for the calculation of the apparent proper motion). The model is based on the standard relativistic-hydrodynamical theory of a relativistic shock that arises from a relativistic explosion into a static, cold external medium (e.g. Meszaros & Rees 1993; Piran et al. 1993; Sari & Piran 1995; Kobayashi et al. 1999; Kobayashi & Zhang 2003; Yi et al. 2013) and is described in detail in Appendix D. We note that the model does not include the possible sideways expansion of the shock.

The free parameters of the model are the energy-to-density ratio E/A, the duration T, the initial Lorentz factor Γ₀, the jet half-opening angle θ_j, the external medium density profile slope k and the viewing angle θ_v. Hereafter we fix T = T₉₀/(1 + z) = 251 s, where T₉₀ refers to the time encompassing the 5% to 95% percentile of the total photon counts as seen in the observer’s reference frame, and we assume Γ₀ = 10³ based on the lower limits from Lesage et al. (2023). Moreover, we consider only two values of the external medium density profile slope, which are k = 0 (homogeneous external medium) and k = 2 (wind-like external medium).

The model predicts the time evolution of the projected angular diameter of the forward and reverse shock images (Eq. D.7). On the other hand, our estimated source sizes are obtained by fitting a circular Gaussian model to the visibility data. The ratio of the two sizes ξ = s_m/s (where s_m is the angular diameter predicted by the model, and s is the FWHM of the circular Gaussian) depends on the lowest-order terms of the MacLaurin expansion in UV radius of the Fourier transform of the surface brightness distribution of the source (Thompson et al. 2017). Taking the surface brightness from Blandford & McKee (1976) as computed in Granot (2008) as reference, we expect 1.2 ≲ ξ ≲ 1.8. This range accommodates values previously considered in the literature, such as the value ξ = 1.4 used in Pihlström et al. (2007) and ξ = 1.3 used in Salafia et al. (2022). Since our model does not predict the surface brightness distribution, we include ξ in our model as a nuisance parameter, with a uniform prior in the range [1.2, 1.8].

4. Results

4.1. Source size expansion

Figure 1 shows the source size constraints from Table 1 in the form of a ‘violin plot’, with the width of each shaded region being proportional to the posterior probability density of the FWHM, horizontally centred at the time of the observation. Additionally, we show the median and 68% symmetric credible interval on the FWHM by means of an error bar for each observation, except for cases where the posterior probability density does not show a clear peak, for which we show instead the 3σ upper limit with a downward-pointing triangle. In order to quantify the source size evolution from these observations, we fit a simple phenomenological power law evolution model, s_m(t_obs)∝t_obs^a, to these size measurements, through the method outlined in Sect. 3.2. The resulting posterior probability density of the power law slope is shown in the inset of Figure 1. The median and symmetric 68% credible interval is $a=0.{69}_{-0.14}^{+0.13}$. We found that more than 99.99% of the posterior probability (> 4 σ-equivalent) is located at a > 0. Therefore, our observations strongly support the expansion of the source. In the main panel of Figure 1, we show with a black dashed line the median of the posterior predictive distribution, that is, the probability distribution of s_m(t_obs) at each fixed t_obs, as derived from the fit. The dotted lines encompass the 68% symmetric credible interval of the same distribution, filled with a grey shade. A lighter grey shading shows the 95% symmetric credible interval. We note that the size measurements in our 15 GHz VLBA epochs at 44 and 205 days are in mild tension with the EVN measurements at similar times. To explore the possibility of a frequency-dependent size, we repeated the power law size evolution model fit considering only observations performed with the EVN or VLBA. Fig. 2 shows the resulting size evolution as fitted to EVN observations at 4.9 GHz and 8.3 GHz (upper panel) or VLBA observations at 15 GHz (lower panel). The plots are similar to Figure 1, except that the epochs not considered in the fit are shown with a light grey shading. The constraint on a from these fits results in medians and symmetric credible intervals of $a=0.{79}_{-0.23}^{+0.19}$ (4.9–8.3 GHz) and $a=0.{98}_{-0.38}^{+0.36}$ (15 GHz), in agreement with each other. On the other hand, the normalisations of the EVN and VLBA power laws differ at the ∼2σ level, as can be evinced from the two-dimensional posterior probabilities shown in Fig. 3.

Fig. 2. Size evolution considering only observations from a single array (upper panel: EVN; lower panel: VLBA). Each panel is similar to Figure 1, except that the epochs not considered in the fit are shown with light grey shading for clarity.

Fig. 3. Comparison of the slope and normalisation of the size evolution as probed by the EVN and the VLBA. The contours in the plot contain 68% (darker contours) and 95% (lighter contours) of the posterior probability on the two parameters (slope and size at a reference time of 100 days) of a single power law fitted to the EVN (blue) or VLBA (orange) size evolution.

In order to exclude the possibility that our results with the EVN are driven by systematic effects, we carried out a series of tests including the check source J1905+1943. We present the results of our tests in Appendix C. The results of these tests indicate that the observed evolution is not driven by systematic errors in the calibrations.

4.2. Apparent proper motion

VLBI observations can constrain the apparent proper motion of the centroid of the emission and, therefore, the jet viewing angle. The source position at each VLBA epoch is displayed in Fig. 4: our results do not show any significant apparent proper motion between 44 and 262 days post-burst, but our statistical errors can accommodate a displacement of up to about 0.6 mas (at the one-σ level) over that period. As shown in Appendix D.2, such an upper limit does not constrain strongly θ_v, which can still be several degrees off the edge of the jet, unless the energy-to-density ratio of the explosion is very large. Still, a number of studies including LHAASO Collaboration (2023) and O’Connor et al. (2023) have used their data to justify a very small θ_v for GRB 221009A, indicating that we are viewing the jet close to on-axis. The lack of significant proper motion observed during our VLBI campaign is fully consistent with such on-axis scenario.

Fig. 4. Source position in our VLBA observations. For each epoch, we show the statistical error bar centered on the median position from the circular Gaussian source fit, with bars spanning the symmetric 68% credible interval on the source position in each direction.

We note that the EVN campaign was not used for such study because of the change in phase reference source between the second and third epoch. While this change was motivated by the discovery of a closer phase calibrator (and hence a more efficient observing strategy), the different systematics and the lack of a reliable a priori position of the new calibrator prevent a reliable astrometric characterisation.

5. Discussion

5.1. Slope of the size evolution with time

The size evolution power law slope $a=0.{69}_{-0.14}^{+0.13}$ we derived is compatible with the expected slopes for a spherical Blandford & McKee (1976) blastwave expanding into a homogeneous medium, α = 5/8 = 0.625, or a wind-like medium, α = 3/4 = 0.75. Moreover, the evolution of the projected physical size is quite similar to that of GRB 030329, the only other burst to date with a measured expansion rate (Taylor et al. 2004, Fig. 5).

Fig. 5. Comparison of size measurements from VLBI observations of GRB 030329 (red circles; Taylor et al. 2004; Pihlström et al. 2007) and GRB 221009A (dark blue squares, this work). Triangles represent upper limits.

On the other hand, at the time of our observations, the anisotropy of the shock due to the finite opening angle of the jet should be observable, both in terms of a steepening (a ‘jet break’, Rhoads 1997) in the light curves, and in terms of a flattening in the evolution of the projected size (Granot et al. 2005). The presence of a jet break in the very high energy afterglow light curve at < 1000 seconds post-trigger has been discussed by LHAASO Collaboration (2023). A jet-break was also suggested by Levan et al. (2023) at ≲2600 seconds post-trigger, using optical to mid-IR data. The expected post-jet-break size evolution slope in the case of a homogeneous external medium is a = 1/4, in absence of an efficient sideways expansion of the shock (Granot et al. 2005). Such a shallow slope is in tension with the observed one at the ∼3σ level. Conversely, the expected average slope is steeper and lies in the range ⟨a⟩∼0.6 − 0.8 if the shock expands sideways (Granot et al. 2005), which is compatible with the observed one. Therefore, in the homogeneous external medium scenario, our observations indicate that either the jet break has not happened yet, or that the shock is expanding sideways. Given the very large isotropic equivalent energy in the gamma-rays, the ‘late jet break’ scenario would pose very demanding requirements on the total energy (see e.g. O’Connor et al. 2023). On the other hand, numerical simulations of external shocks arising from relativistic jets and analytical arguments seem to indicate that the sideways expansion is inefficient, unless the initial opening angle is very narrow (van Eerten et al. 2010; De Colle et al. 2012; Granot & Piran 2012). These difficulties could be alleviated if the jet features a structure consisting of a narrow ‘core’ surrounded by ‘wings’ where the kinetic energy per unit solid angle decreases slowly, as suggested by O’Connor et al. (2023) and Gill & Granot (2023). This profile would steepen the evolution of the observed size, making it more similar to the spherical case, but with a reduced energy requirement with respect to a wide jet with a uniform angular energy profile.

In the wind medium case, the expected post-jet-break size evolution slope is a = 1/2 in absence of sideways expansion. Therefore, in such a scenario, the observed evolution does not indicate the need for sideways expansion nor for a non-uniform structure within the opening angle.

5.2. Possible frequency-dependent size

As discussed in section 4.1, our data suggest the presence of a frequency-dependent size evolution. We explore here a possible avenue to interpret this behaviour. The radio afterglow of this GRB cannot be explained by a simple FS propagating either into a wind-like or a homogeneous environment (Ren et al. 2023; Sato et al. 2023; Ren et al. 2024; Zheng et al. 2024). Using a data set encompassing observations from the GeV to the radio domain, Laskar et al. (2023) showed that the standard afterglow model struggles at explaining the radio emission both with a FS and a RS of a conical jet propagating through a wind-like environment, leading them to invoke an additional component whose temporal evolution does not follow the standard prescriptions. Such a component could be a RS, dominating the emission at the lower frequencies (≲10 GHz) up to t_obs ≲ 100 d (see the modelling of O’Connor et al. 2023; Gill & Granot 2023).

Stimulated by these studies, we explored a scenario where the emission we observed is a superposition of a FS and a RS. We fitted the model described in Sect. 3.3 to the observed size evolution, leaving E/A and θ_j as our free parameters, and additionally including the ξ nuisance parameter (see section 3.3). The external medium power law index was fixed to k = 0 or k = 2. Based on the apparent possible frequency-dependent behaviour, we assumed the higher-frequency observations (15 GHz) to be dominated by the FS, while the lower-frequency ones (4.9–8.3 GHz) to be dominated by the RS. For each external density profile, the posterior probability density of (E/A, θ_j, ξ) was derived through the Bayesian approach described in Sect. 3.2, with uniform-in-log priors on E/A in the range [10⁵², 10⁵⁸] erg cm³ and on θ_j in the range [0.5, 30] deg, and a uniform prior on ξ in the range [1.2, 1.8], as explained in section 3.3. Figure 6 shows the FS and RS model size evolution fitted to the observations. Figure 7 shows corner plots of the posterior probability densities, marginalized over ξ.

Fig. 6. Model size evolution in a FS plus RS scenario. The violins and the error bars show the source size evolution as inferred by our observations, in the same way as in Figure 1. The gray dotted line and orange dashed lines show the medians of the posterior predictive distributions of the FS and RS size, respectively, as obtained by fitting the physical model described in Appendix D to the sizes shown in the figure, assuming a homogeneous (top panel) or wind-like (bottom panel) external medium, and assuming the FS to dominate at 15 GHz and the RS to dominate at 4.9 and 8.3 GHz. The shaded bands around these lines show the 68% credible interval of the posterior predictive distribution.

Fig. 7. Corner plot of the posterior probability density of the physical model parameters in the homogeneous (top panel) or wind-like (bottom panel) external medium case and assuming the FS to dominate at 15 GHz and the RS to dominate at 4.9–8.3 GHz. In each panel, the red histograms show the marginalised posterior probability densities, with black solid vertical lines showing the median and dashed lines showing the 68% credible interval or, if the latter extends to the lower (upper) extremum of the prior range, the 95% upper (lower) limit (values reported on top of the panels). The filled contours show the smallest regions containing 68% and 95% of the two-dimensional posterior probability, with the black squares showing the position of the median. The grey lines show contours of constant total jet energy assuming A = 1 cm−3. Each contour is labeled with the base-10 logarithm of the corresponding total jet energy.

The model can reasonably reproduce the data in both the homogeneous (k = 0) and wind-like (k = 2) external medium cases. Due to the relatively steep observed size evolution, and given that our model does not include sideways expansion, the required jet half-opening angle is large (θ_j > 13 deg at the 95% credible level for k = 0; θ_j > 16 deg for k = 2). Assuming the standard reference external density A = 1 cm⁻³, this in turn pushes the total jet energy E_jet = E(1 − cos θ_j) to very large values, as demonstrated by the grey contours in Fig. 7. The energy requirement can be reduced if the external density is much lower than the reference value, A ≪ 1 cm⁻³, or if the jet features an angular structure such that the energy per unit solid angle decreases away from the jet axis, as discussed in the previous section. In this case, the decrease must be shallow (E ∝ θ^−a with a < 2, where θ is the angle from the jet axis), otherwise both the emission and the size evolution would simply reflect those of a uniform jet (O’Connor et al. 2023), leading again to the same difficulties with the jet opening angle and total energy.

6. Conclusions

In this paper, we presented VLBI observations of the brightest γ-ray burst ever observed, GRB 221009A. The high angular resolution provided by the EVN and the VLBA allowed us to constrain the size and the expansion of the blast wave produced by the GRB ejecta for the second time ever. The expansion rate is consistent with the expectation for a spherical Blandford & McKee (1976) blast wave (i.e., an ultra-relativistic FS) propagating into a homogeneous or wind-like medium. This could be taken as an indication that the shock is anisotropic only on angular scales larger than those probed by our observations (i.e. the jet break has not happened yet). This in turn points to an extremely large total (collimation-corrected) energy in the jet, especially if the external medium features a homogeneous density. The demanding energy requirement could be alleviated if the external density were much lower than usually assumed, or if the jet energy per unit solid angle decreased slowly with the angle from the jet axis, as proposed by O’Connor et al. (2023) and Gill & Granot (2023). Alternatively, the shock could be undergoing sideways expansion during the time covered by our observations.

Additionally, our observations suggest a frequency-dependent size evolution, with the VLBA observations at 15 GHz showing a somewhat larger size at 40 days after the explosion, and a faster increase afterwards, with respect to the EVN observations at 5 and 8 GHz. This could be due to the emission being dominated by the reverse shock at the lower frequencies, and by the forward shock at the higher frequencies. Our work highlights the crucial role played by multi-wavelength VLBI monitoring of transient events both at early and late times, and in providing a vital insight into the physics of such events.

Data availability

A copy of the reduced visibilities is available at the CDS via anonymous ftp to cdsarc.cds.unistra.fr (130.79.128.5) or via https://cdsarc.cds.unistra.fr/viz-bin/cat/J/A+A/690/A74

Acknowledgments

The European VLBI Network is a joint facility of independent European, African, Asian, and North American radio astronomy institutes. Scientific results from data presented in this publication are derived from the following EVN project code: RG013. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. This work made use of the Swinburne University of Technology software correlator, developed as part of the Australian Major National Research Facilities Programme and operated under licence. SG would like to thank Z. Paragi and the staff of JIVE for their help and support during his visiting period in Dwingeloo. We would like to thank the directors and staff of all the EVN telescopes for approving, executing, and processing our out-of-session ToO observations. The research leading to these results has received funding from the European Union’s Horizon 2020 Research and Innovation Programme under grant agreement No. 101004719 (OPTICON RadioNet Pilot). The research leading to these results has received funding from the European Union’s Horizon 2020 Programme under the AHEAD2020 project (grant agreement n. 871158). This work has been funded by the European Union-Next Generation EU, PRIN 2022 RFF M4C21.1 (202298J7KT – PEACE). OS acknowledges funding from the Istituto Nazionale di Astrofisica, project number 1.05.23.04.04. BM acknowledges financial support from the State Agency for Research of the Spanish Ministry of Science and Innovation under grant PID2019-105510GB-C31/AEI/10.13039/501100011033 and through the Unit of Excellence María de Maeztu 2020–2023 award to the Institute of Cosmos Sciences (CEX2019- 000918-M). MPT acknowledges financial support through grants CEX2021-001131-S and PID2020-117404GB-C21 funded by the Spanish MCIN/AEI/ 10.13039/501100011033.

References

Appendix A: EVN observation strategy and data reduction

In this Appendix we provide detailed information on the observation strategy and the data reduction of the EVN observations. As explained in Section 2, the structure of the observations followed a typical phase-referencing experiment. Three compact, extragalactic radio sources J1800+3848, J1925+2106 and J0121+0422 were used as fringe finders and bandpass calibrators in the campaign. The target scans, lasting approximately 4.5 and 2.5 minutes at 5 and 8 GHz, respectively, were interleaved with 1.5 minute scans of the phase calibrator. In the first two observations, namely 40 and 43 days post-burst (RG013 B and C), the radio source J1905+1943 was used as a phase calibrator and the Very Large Array Sky Survey (VLASS) compact radio source J191142+1952 was included for testing its suitability as a closer phase reference source (d = 0.33° from the GRB position). Given the success of that test, J191142+1952 was then adopted as phase calibrator in the last three epochs, i.e. from 117 to 261 days post-burst (RG013 D, E and F). In order to inspect the consistency of the calibration procedure, one or multiple compact radio sources were observed approximately every 30 minutes. The phase and amplitude solutions derived from the calibrators were applied to these check sources to verify the quality of the calibration.

The calibration was performed using AIPS (Greisen 2003), following the standard procedure for EVN phase-referenced observations. The amplitude calibration, which accounts for the bandpass response, the antenna gain curves and the system temperatures, was performed using the results from the EVN pipeline. Procedures vlbatecr and vlbampcl were used to correct for the dispersive delay and to calculate the manual single band delay on the fringe finder, respectively. Subsequently, we carried out the global fringe fitting on the phase calibrator with the task fring. Solutions were interpolated and applied to the phase calibrator itself, the check sources and the target with the task clcal. After this point, the calibration procedure differs according to the epoch.

For the first epoch (RG013 B), we carried out the fringe fitting on J1905+1943 using a model of the source derived by a concatenation (in CASA, McMullin et al. 2007) and self-calibration (in Difmap) of all the visibilities on the source obtained across the two epochs at 8 GHz. This approach is warranted by the stability of the structure of extragalactic sources on the duration of the campaign, and improves the phase, delay, and rate calibration by accounting for the possibile structure of the phase-reference source. We then interpolated the solutions and applied the results to J1905+1943 itself, the check source (J1923+2010) and GRB 221009A. Moreover, for RG013 B and C, we tried two rounds of self-calibration on J1905+1943, first in phase-only, with a solution interval of 2 minutes, and then in amplitude and phase with a solution interval of 60 minutes in AIPS. Since the self-calibration on J1905+1943 did not result in a significant improvement in the final S/N of the image of the GRB, we used the image generated with the calibration performed prior to the self-calibration.

In the last three epochs, we employed a different phase calibrator, J191142+1952, motivated by the significantly smaller separation from the source (0.33° vs 1.75°). If the phase calibrator is closer to the target, any possible decorrelation of the phase solutions is significantly reduced. However, the position of J191142+1952 was constrained only with a precision of the order of an arcsecond: for VLBI observations, this means that the coordinates of the centre of the source were not aligned with the phase centre of the observation. If one does not correct for the uncertainty in the position, the phase solutions of the global fringe fitting on the phase calibrator will contain a systematic error and the apparent coordinates of the centre of the sources to which these solutions are applied will be incorrect. To avoid this, we started from the fourth epoch, made at 8 GHz, which provides higher angular resolution and therefore a more precise position of the calibrator. We applied the solutions of the first global fringe fitting on J191142+1952 to the check sources, J1905+1943 and J1923+2010, we produced an image of each of them and we compared the apparent coordinates of each source with the actual position, known with an uncertainty of the order of mas. Since J1905+1943 and J1923+2010 appear to be aligned in the sky, with J191142+1952 placed in between, we derived the real coordinates with a 1D interpolation at the position of J191142+1952 of the offset observed for the check sources. We then re-calculated the visibilities of J191142+1952 by fixing the phase centre with the fixvis task in CASA, inserting the new coordinates. We then repeated the entire calibration process iteratively, until the apparent and the real sky coordinates of the check sources were consistent within the resolution of the observation. We corrected the third and fifth epoch using the position derived at 8.3 GHz.

Subsequently, we produced a model of J191142+1952 in Difmap and we used the model as input to perform the global fringe fitting of the phase calibrator in AIPS, in order to take into account any possible structure of J191142+1952 and correct for it. We interpolated the solutions and we applied them to J191142+1952, J1905+1943, J1923+2010 and GRB 221009A. Lastly, we performed a round of amplitude and phase self-calibration of J191142+1952 in AIPS, using a solution interval of 2 minutes and we applied the interpolated solutions to J191142+1952, J1905+1943, J1923+2010 and GRB 221009A. After each of the aforementioned steps of the procedure, the derived solutions were inspected and bad data were properly flagged.

Images of the target and the check sources were produced using Difmap. For the RG013 C dataset, both the noise pattern in the dirty image and the modelfit residuals indicated the presence of some residual phase calibration errors. Of all our datasets, this is the epoch with the largest flux density (S = 1.76 mJy, as estimated interpolating the VLA measurments reported by Laskar et al. 2023). Therefore, considering the presence of the sensitive Effelsberg and Tianma65 radio telescopes and the large bandwidth, we carried out a phase-only self calibration with a solution time interval of one minute and combining all the frequency sub-bands. After carefully checking the behaviour of the solutions, we concluded that the process improved the data quality without introducing any artifact in the data. We carried out the subsequent steps of analysis as for the other epochs.

Appendix B: Circular Gaussian fits to source visibilities

In this appendix we provide more detailed information about our circular Gaussian fits to source visibilities from our VLBI observations. Our observations and calibration procedure yielded a set of N complex visibilities {𝒱_i}_i = 1^N, which measure the two-dimensional Fourier transform of the image at positions (u_i, v_i) in the (u, v) plane, with associated uncertainties σ_i = w_i^−1/2, where w_i are the ‘data weights’ determined by the calibration procedure. Each visibility resulted from interferometry between the pair of antennae (‘baseline’) B_i = (b_1, i, b_2, i), where b_1, i and b_2, i are indices that identify the two antennae, which satisfy 1 ≤ b_1, i < N_A and b_1, i < b_2, i ≤ N_A, were N_A is the total number of antennae in the VLBI network. The total number of distinct baselines is N_B = N_A(N_A − 1)/2.

We assumed a Gaussian likelihood ℒ for these visibilities,

$$\begin{array}{c}\hfill ln\mathcal{L}\left(\{{\mathcal{V}}_i{\}}_{i=1}^N|\mathit{x}\right)=-\frac{1}{2}{\displaystyle \sum _{i=0}^N}\frac{{({\mathcal{V}}_{\mathrm{m}}(u_i,v_i,\mathit{x})-{\mathcal{V}}_i)}^2}{{\sigma }_i^2},\end{array}$$(B.1)

where x is the model parameter vector and 𝒱_m(u_i, v_i, x) is the visibility predicted by the model. We adopted a circular Gaussian surface brightness model, with amplitude correction factors introduced to account for potential antenna-based systematic mis-calibrations (similar to Natarajan et al. 2017), namely

$$\begin{array}{c}\hfill {\mathcal{V}}_{\mathrm{m}}(u_i,v_i,\mathit{x})=G_{b_{1,i}}^{-1}{G}_{b_{2,i}}^{-1}{F}_{\nu }{e}^{-2{\pi }^2{\left(\frac{s}{\sqrt{8ln2}}\right)}^2(u_i^2+v_i^2)-2\pi j(u_i\rho +v_i\delta )},\end{array}$$(B.2)

where $j=\sqrt{-1}$, F_ν is the total flux density, s is the full width at half maximum of the circular Gaussian, ρ and δ are the spherical offsets of the source with respect to the phase centre, and G_b1, i and G_b2, i are dimensionless factors that encode the amplitude mis-calibrations of antennae b_1, i and b_2, i (Natarajan et al. 2017). These are defined as follows: G_i = 1 means that there is no systematic error in the amplitude calibration of antenna i; conversely, G_i = 0.9 (resp. G_i = 1.1) means that the amplitude of that particular antenna is systematically overestimated (resp. underestimated) by 10%.

The parameter vector of the model is therefore x = (F_ν, s, ρ, δ, G₁, ..., G_NA). By Bayes’ theorem, we defined the posterior probability on x, given our data $\{{\mathcal{V}}_i{\}}_{i=1}^N$, as

$$\begin{array}{c}\hfill P\left(\mathit{x}|\{{\mathcal{V}}_i{\}}_{i=1}^N\right)\propto \pi (\mathit{x})\mathcal{L}\left(\{{\mathcal{V}}_i{\}}_{i=1}^N|\mathit{x}\right),\end{array}$$(B.3)

where π(x) is the prior probability on the parameters. We decomposed the latter into a product of independent one-dimensional priors on all parameters. We adopted simple independent uniform priors on the source parameters, with the due bounds F_ν > 0 and s > 0. Where necessary, in order to prevent the fitting procedure from picking up some noise peak instead of the actual source, we restricted the position (ρ, δ) to within a small angular distance Δ_pos ∼ 1 mas from the peak (ρ₀, δ₀) of the cleaned map constructed with Difmap. All, but one, of the priors on the mis-calibration parameters G_i were set equal to Gaussians centered at 1 with sigma equal to 0.1, hence admitting typical systematic errors of 10%. The prior on the remaining G_i parameter was set to a delta function centered at 1, in order to break a degeneracy with the total flux density (see, e.g., Natarajan et al. 2017). The index i of such fixed calibration parameter corresponds to the reference antenna in the network, which we chose to be Effelsberg for the EVN and Fort Davis for the VLBA (with the exception of epoch BA160C1, where we used Los Alamos).

For each epoch, we sampled the posterior probability using the EMCEE (Foreman-Mackey et al. 2013) python package. We initialised EMCEE with the initial guess x = (I_ν, pk, 10⁻² mas, ρ₀, δ₀, 1, ..., 1), where I_ν, pk is the peak surface brightness (expressed in Jy/beam) in the cleaned map, corresponding to an unresolved circular Gaussian source at the position of the peak of the cleaned map and with a flux density that yields the observed peak surface brightness. We then ran 5000 iterations of the MCMC (depending on the epoch) with 32 walkers, producing ∼10⁵ samples of the posterior probability density for each epoch. The results were constructed after discarding the initial 30% of these samples in each chain as burn-in.

Figure B.1 shows corner plots that illustrate the properties of the posterior probability density of (F_ν, s) for our EVN 4.9 GHz epochs. Figures B.2 and B.3 show the corresponding corner plots for our EVN 8.3 GHz epochs and for our VLBA 15 GHz epochs, respectively.

Fig. B.1. Posterior probability distribution on the source size and flux density in our EVN 4.9 GHz epochs at T − T0 = 43 (RG1013 C, left), 117 (RG013 D, centre) and 261 (RG013 F, right) days. In each corner plot, the top-left and bottom-right sub-panels show histograms of the posterior samples of Fν and FWHM, with the vertical solid lines showing the median and the vertical dashed lines bracketing the 68% credible interval or, if the latter extends to 0, the 95% credible upper limit. The bottom-left sub-panel of each corner plot shows the smallest contours containing 68% and 95% of the posterior probability.

Fig. B.2. Similar to Figure B.1, but for our EVN 8.3 GHz epochs at T − T0 = 40 days (EVN B, left) and 118 days (EVN E, right).

Fig. B.3. Similar to Figure B.1, but for our VLBA 15 GHz epochs at T − T0 = 44 (VLBA B, upper left), 114 (VLBA C, upper right), 205 (VLBA C1, lower left) and 262 (VLBA D, lower right) days.

Appendix C: Tests on the evolution of the flux density and the size

In this Appendix we present tests on the EVN observational results that we carried out in order to exclude the possibility that the measured evolution of the GRB size is a result of systematic effects. These tests include the check source J1905+1943. Unfortunately, due to the sparse (u, v) plane coverage and the large separation, J1923+2010 could not be used to get meaningful constraints. First, the measured GRB afterglow size as a function of the area of the synthesised beam are presented in Fig. C.1. These quantities are clearly not correlated, hence we can exclude the possibility that the observed expansion of the GRB is driven by a systematic change in the width of the synthesised beam. In Fig. C2, the flux density (left panel) and the size (right panel) of GRB 221009A and the check source J1905+1943 are compared. The decrease in the GRB 221009A flux density is not accompanied by a variation of the J1905+1943 flux density, as expected. Regarding the size measurements, while the FWHM of J1905+1943 is constant throughout the campaign, the size of GRB 221009A is clearly increasing in time. Therefore, the variations observed for the GRB afterglow cannot be ascribed to a systematic effect due to an imprecise calibration. Concerning the VLBA, no test was performed because of the lack of close enough check sources.

Fig. C.1. Estimates of the FWHM of the GRB as a function of the total area of the synthesised beam for the EVN observations.

Appendix D: Model of the projected size and proper motion of the forward and reverse shock

D.1. Dynamics and size evolution

In the following, we describe an approximate analytical model of the dynamics of the forward and reverse shocks, based on calculations similar to those of Kobayashi et al. (1999), Kobayashi & Zhang (2003), Yi et al. (2013). The aim is to extend approaches such as those described in Oren et al. (2004) and Granot et al. (2005) by including the reverse shock, which was not considered there. We assume a cold external medium with a power law density profile ρ = Am_p(R/R_⋆)^−k, where R is the radial distance from the progenitor, m_p is the proton mass and A is the number density at a reference radius R_⋆ = 5.5 × 10¹⁷ cm. With this definition, A plays the role of either the homogeneous interstellar medium (ISM) number density, A ≡ n if k = 0, or that of the wind density parameter, A ≡ A^⋆ if k = 2. We assume a simplified description of the jet as a cold, kinetic-energy-dominated shell with uniform initial bulk Lorentz factor Γ₀ = 10³ Γ_0, 3 ≫ 1 and constant isotropic-equivalent kinetic luminosity L = E/T, where E = 10⁵⁵ E₅₅ erg is the isotropic-equivalent jet energy and T = T₉₀/(1 + z)≈251 T_2.4 s (where T_2.4 = T/(10^2.4 s)) is the lifetime of the central engine. The Sedov length associated with this shell is ℓ_S = [(3 − k)E/(4πAR_⋆^km_pc²)]^{1/(3 − k)}. As this shell expands into the external medium at relativistic speed, a FS arises, which sweeps the external medium moving with a Lorentz factor ${\mathrm{\Gamma }}_{\mathrm{FS},0}\sim \sqrt{2}{\mathrm{\Gamma }}_0$. The shocked external medium resides in the region contained between the FS and the contact discontinuity (CD) that separates it from the jet material. Since this implies some deceleration of the jet material behind the CD as well, as soon as the ram pressure of such material overcomes the pressure in the jet (formally already at R = 0 given our assumption of a cold jet), a RS also arises, which separates shocked from cold un-perturbed jet material. Let us indicate with numbers from 1 to 4 the un-perturbed external medium, shocked external medium, shocked jet and un-perturbed jet respectively, as usual. The RS is initially non-relativistic (i.e. the relative speed of regions 3 and 4 is β₃₄ ≪ 1), but it can become relativistic before the RS crosses the whole jet if the condition (Sari & Piran 1995; Kobayashi et al. 1999; Chevalier & Li 2000; Kobayashi & Zhang 2003; Zou et al. 2005; Yi et al. 2013)

$$\begin{array}{c}\hfill \frac{E}{A}<\frac{4\pi m_{\mathrm{p}}{c}^2{R}_{\star }^k}{3-k}{\left(cT\right)}^{3-k}{\mathrm{\Gamma }}_0^{2(4-k)}\approx \{\begin{array}{cc}2.7\times {10}^{60}{\mathrm{\Gamma }}_{0,3}^8{T}_{2.4}^3\mathrm{erg}{\mathrm{cm}}^3& \hfill k=0\\ 4.3\times {10}^{58}{\mathrm{\Gamma }}_{0,3}^3{T}_{2.4}\mathrm{e}\mathrm{r}\mathrm{g}\mathrm{c}{\mathrm{m}}^3& \hfill k=2\end{array}\end{array}$$(D.1)

is satisfied, in which case the jet deceleration is said to be in the ‘thick shell regime’. In the following we describe the dynamics in such regime, and we defer to later the treatment of the opposite, ‘thin shell’ regime. For the homogeneous ISM case, k = 0, the RS transitions from Newtonian to relativistic at a radius $R_{\mathrm{N}}\sim {\ell }_{\mathrm{S}}^{3/2}{[{(12cT)}^{1/2}{\mathrm{\Gamma }}_0^2]}^{-1}\sim 4.2\times {10}^{15}{E}_{55}^{1/2}{A}^{-1/2}{\mathrm{\Gamma }}_{0,3}^{-2}{T}_{2.4}^{-1/2}\mathrm{cm}$, while in the wind case, k = 2, the RS is always relativistic as long as condition D.1 holds. As regions 2 and 3 decelerate due to the increasing amount of swept external medium mass, at some point the RS crosses the whole jet, at a radius

$$\begin{array}{c}\hfill R_{\otimes }={(4(3-k)cT)}^{1/(4-k)}{\ell }_{\mathrm{S}}^{(3-k)/(4-k)}.\end{array}$$(D.2)

Before R_⊗, regions 2 and 3 effectively expand at the same Lorentz factor Γ, whose evolution can be described approximately as

$$\begin{array}{c}\hfill \mathrm{\Gamma }(R)\sim \{\begin{array}{cc}{\mathrm{\Gamma }}_0\hfill & \hfill R\le R_{\mathrm{N}}\\ \frac{{\ell }_{\mathrm{S}}^{(3-k)/4}}{{[4(3-k)cT]}^{1/4}}{R}^{-(2-k)/4}\hfill & \hfill R_{\mathrm{N}}<R\le R_{\otimes }\end{array}.\end{array}$$(D.3)

The Lorentz factor of region 3 at the end of the RS crossing is therefore Γ_⊗ = [ℓ_S/(4(3−k)cT)]^{(3 − k)/[2(4 − k)]}. At radii larger than R_⊗, the Lorentz factor of region 2 follows the Blandford & McKee (1976) relativistic blastwave evolution, Γ₂ ∼ (R/ℓ_S)^{−(3 − k)/2}. This holds as long as the lateral expansion of the shocked material in region 2 is negligible: numerical simulations and analytical arguments (e.g. Kumar & Granot 2003; van Eerten et al. 2010; De Colle et al. 2012; Lyutikov 2012; Granot & Piran 2012) show that such expansion has a very limited impact on the dynamics until region 2 becomes mildly relativistic, which justifies such assumption. In the homogeneous ISM case, k = 0, the subsequent evolution of Γ₃ has been historically described phenomenologically (Kobayashi & Sari 2000) as Γ₃ = Γ_⊗(R/R_⊗)^−g, with g being typically fixed at around g ∼ (7 − 2k)/2 in the case of a non-relativistic RS, or at g = (3 − k)/2 (Kobayashi et al. 1999; Yi et al. 2013) in the case of a relativistic RS (when condition D.1 holds), based on insights from the numerical simulations described in Kobayashi et al. (1999) and Kobayashi & Sari (2000). Physically, the different evolution is likely related to the conversion of internal to kinetic energy in region 3 as it expands, which allows it to remain ‘attached’ to region 2 as long as its temperature is relativistic. For historical reasons, in the case of a wind environment, k = 2, the evolution in this phase has been always assumed to track that of the tail of the FS Blandford-McKee solution (i.e. g = −3/2), despite the lack of numerical simulations to compare to. We argue here that generally, as long as region 3 contains more internal than kinetic energy, its acceleration will keep it ‘attached’ to region 2, and they will evolve following the Blandford & McKee (1976) solution. As the internal energy conversion terminates, region 3 must eventually ‘detach’ and expand backwards (as seen from the CD) into a rarefaction wave, and thus the evolution of Γ₃ with radius must steepen. In order to estimate the radius R_det at which regions 2 and 3 detach, we need to know the evolution of the internal energy in region 3, E_int, 3 (as measured in the comoving frame of region 3). From the first equation of thermodynamics, $dln{E}_{\mathrm{int},3}=(1-\widehat{\gamma })dln{V}_3\prime$, where $\widehat{\gamma }$ is the adiabatic index and V₃′ is the comoving volume of region 3. We assume V₃′∝R³/Γ₃. Right after the shock crossing regions 2 and 3 still move together, hence we can assume Γ₃ ∝ R^{−(3 − k)/2}, which leads to $E_{\mathrm{int},3}=E_{\mathrm{int},3,\otimes }{(R/R_{\otimes })}^{(1-\widehat{\gamma })(9-k)/2}$. Taking the internal energy at the end of RS crossing to be E_{int, 3, ⊗} ∼ (Γ_34, ⊗ − 1)m₃c² ∼ (Γ_⊗/Γ₀ + Γ₀/Γ_⊗)m₃c²/2 (where Γ_34, ⊗ is the relative Lorentz factor of regions 3 and 4 at the RS crossing radius, and m₃ is the jet rest mass), we finally conclude that the effective dimensionless temperature in region 3 evolves as

$$\begin{array}{c}\hfill {\mathrm{\Theta }}_3\equiv E_{\mathrm{int},3}/m_3{c}^2\sim [\frac{1}{2}(\frac{{\mathrm{\Gamma }}_{\otimes }}{{\mathrm{\Gamma }}_0}+\frac{{\mathrm{\Gamma }}_0}{{\mathrm{\Gamma }}_{\otimes }})-1]{(R/R_{\otimes })}^{(1-\widehat{\gamma })(9-k)/2}.\end{array}$$(D.4)

Assuming $\widehat{\gamma }=4/3$ since the RS is relativistic, we finally obtain the detachment radius from the condition Θ₃(R_det) = 1, which yields

$$\begin{array}{c}\hfill R_{\det }=\max \{{[\frac{1}{2}(\frac{{\mathrm{\Gamma }}_{\otimes }}{{\mathrm{\Gamma }}_0}+\frac{{\mathrm{\Gamma }}_0}{{\mathrm{\Gamma }}_{\otimes }})-1]}^{6/(9-k)},1\}{R}_{\otimes },\end{array}$$(D.5)

where the maximum function is introduced to account for cases where Θ₃ < 1 at R_⊗, in which case R_det = R_⊗.

Based on these considerations, we model the evolution of Γ₃ after RS crossing as

$$\begin{array}{c}\hfill {\mathrm{\Gamma }}_3(R)=\{\begin{array}{cc}{\mathrm{\Gamma }}_{\otimes }{\left(\frac{R}{R_{\otimes }}\right)}^{-(3-k)/2}\hfill & \hfill R_{\otimes }\le R<R_{\det }\\ {\mathrm{\Gamma }}_{\otimes }{\left(\frac{R_{\det }}{R_{\otimes }}\right)}^{-(3-k)/2}{\left(\frac{R}{R_{\det }}\right)}^{-(7-2k)/2}\hfill & \hfill R_{\det }\le R\end{array}.\end{array}$$(D.6)

The above relations completely specify the evolution of Γ₂ and Γ₃ with radius as a function of the Sedov length ℓ_S, initial Lorentz factor Γ₀ and jet duration T for a given choice of k in the thick shell regime. The thin shell regime is obtained (Kobayashi et al. 1999) by setting all transition radii equal to the ‘deceleration’ radius, R_N = R_⊗ = R_det = ℓ_S/Γ₀^{2/(3 − k)}. The relation between the radius R and the observer time for region i ∈ {2, 3} can be obtained by noting that most of the emission that the observer receives comes from material moving at an angle ∼Γ_i⁻¹ from the line of sight, for which the arrival time is t_obs/(1 + z) = t(R)−Rβ_i(R)/c, where ${\beta }_i=\sqrt{1-{\mathrm{\Gamma }}_i^{-2}}$. The progenitor-frame time t as a function of the radius can be obtained by integrating t(R) = ∫₀^RdR/(β_ic). By numerically inverting the resulting relation between t_obs and R, we finally obtain R_i(t_obs), and thus Γ₂(t_obs) and Γ₃(t_obs). The projected angular diameter of region i = 2, 3 is then approximately (e.g. Oren et al. 2004)

$$\begin{array}{c}\hfill s_{\mathrm{m},i}(t_{\mathrm{obs}})\sim 2\frac{R_i(t_{\mathrm{obs}})}{d_{\mathrm{A}}}\times \{\begin{array}{cc}{\mathrm{\Gamma }}_i^{-1}(t_{\mathrm{obs}})\hfill & \hfill {\mathrm{\Gamma }}_i(t_{\mathrm{obs}})\ge 1/{\theta }_{\mathrm{j}}\\ {\theta }_{\mathrm{j}}\hfill & \hfill {\mathrm{\Gamma }}_i(t_{\mathrm{obs}})<1/{\theta }_{\mathrm{j}}\end{array},\end{array}$$(D.7)

where θ_j is the jet half-opening angle. We find that the predicted size of the FS from this model matches that of the more refined model of Granot et al. (2005) within 10% for k ∈ {0, 2} in the self-similar deceleration stage.

D.2. Proper motion

Unless the jet is observed perfectly on-axis, some apparent displacement of the projected image is expected. We discuss here the case θ_v > θ_j, which produces the largest such displacement. As long as Γ⁻¹ < (θ_v − θ_j), the observed emission is dominated by the border of the shock closest to the observer. Its apparent displacement Δ can be modelled effectively as that of a point source moving at ∼c at an angle θ_v − θ_j away from the line of sight, so that the displacement increases linearly in time, Δ ∝ t_obs. For (θ_v − θ_j) < Γ⁻¹ < (θ_v + θ_j), the emission is dominated by material moving at ∼1/Γ from the line of sight, hence the displacement evolves as Δ(t_obs)∼R(t_obs)/Γ(t_obs)d_A, with R(t_obs) being the same as in the on-axis case. Eventually, for Γ⁻¹ > (θ_v + θ_j), the emission is dominated by the material at the shock border farthest from the observer, and the displacement is therefore Δ(t_obs)∼R(t_obs)sin(θ_v + θ_j)/d_A. The transition times t_{j, obs, ±} that separate the three regimes described above can be obtained by setting Γ(t_obs) = (θ_v ± θ_j)⁻¹ in the on-axis case, which we do numerically.

The model described here neglects the effects of lateral expansion of the shock and of a non-trivial jet structure outside the ‘core’ of half-opening angle θ_j. The former would generally slow down the evolution, so that the displacement predicted by this model can be considered as an upper limit. The latter would change (generally steepen) the slope of the displacement with time before the time t_{j, obs, −} at which the jet core starts coming into sight, but not thereafter. For the most likely parameters, our observations are at t_obs > t_{j, obs, −}, so that the effects of a jet structure are unimportant for this particular source.

Figure D1 shows the displacement predicted by such a model between 44 d and 262 d, assuming the emission to be dominated by the FS (which produces the largest displacement, and dominates the VLBA data according to our interpretation) for different assumptions on θ_j (varying across columns) and on the external medium power law index (top row: k = 2; bottom row: k = 0), as a function of the off-edge viewing angle θ_v − θ_j and of the energy to density ratio E/A. These predictions show that our upper limit on the observed displacement only excludes off-edge viewing angles between a few degrees and around 11 degrees, combined with large energy to density ratios E/A ≳ 10⁵⁵ erg cm³ for k = 2, or rather extreme E/A ≳ 10⁵⁸ erg cm³ for k = 0. Viewing angles larger than θ_v − θ_j ∼ 11 deg cannot be constrained because in that case the shock is in the ‘point-source at the jet edge’ regime all the way to 262 d, with a rather small apparent transverse velocity. On the other hand, such a large viewing angle would be very unlikely given the huge γ-ray isotropic-equivalent energy of this source.

Fig. D.1. Flux density (left panel) and size (right panel) of GRB 221009A and J1905+1943 in the EVN observations. The average flux density of J1905+1943 (left panel) and the 1:1 correlation (right panel) are shown as grey dashed lines for the sake of comparison.

Fig. D.2. Constraint on the viewing angle from the absence of an observed source apparent displacement in our VLBA observations. In each panel, filled contours show the displacement of the centre of the fitted Gaussian expected between 44 d and 262 d, color coded as shown in the colorbar on the right, as a function of the E/A ratio and of the off-edge viewing angle θv − θj. The red contour shows Δ(262 d)−Δ(44 d) = 0.6 mas, which represents the largest displacement compatible at 1 σ with our observations. The red contour hence contains the excluded parameter region. The upper panel row refers to a wind-like external medium, while the lower row refers to a homogeneous external medium. Each column assumes a different jet half-opening angle, given at the top of the column.

All Tables

All Figures

Fig. 1. Source size as a function of time. The source size constraints obtained as described in Section 3.1 are shown in the form of violin plots of different colours, centred at the observing time of each epoch and of proportional width to the posterior probability density of the FWHM. In addition, we show the median and 68% credible interval with an error bar of the same color, or the 95% credible upper limit with a triangle if the former interval extends to 0. The black dashed line and the two grey shaded areas show respectively the median, 68% credible interval and 95% credible interval of the posterior predictive distribution of the source size evolution obtained from fitting a power law model s ∝ tobsa to the sizes from all the epochs. The inset shows the posterior probability density of the slope a from such fit.

In the text

	Fig. 2. Size evolution considering only observations from a single array (upper panel: EVN; lower panel: VLBA). Each panel is similar to Figure 1, except that the epochs not considered in the fit are shown with light grey shading for clarity.
In the text

	Fig. 3. Comparison of the slope and normalisation of the size evolution as probed by the EVN and the VLBA. The contours in the plot contain 68% (darker contours) and 95% (lighter contours) of the posterior probability on the two parameters (slope and size at a reference time of 100 days) of a single power law fitted to the EVN (blue) or VLBA (orange) size evolution.
In the text

	Fig. 4. Source position in our VLBA observations. For each epoch, we show the statistical error bar centered on the median position from the circular Gaussian source fit, with bars spanning the symmetric 68% credible interval on the source position in each direction.
In the text

	Fig. 5. Comparison of size measurements from VLBI observations of GRB 030329 (red circles; Taylor et al. 2004; Pihlström et al. 2007) and GRB 221009A (dark blue squares, this work). Triangles represent upper limits.
In the text

Fig. 6. Model size evolution in a FS plus RS scenario. The violins and the error bars show the source size evolution as inferred by our observations, in the same way as in Figure 1. The gray dotted line and orange dashed lines show the medians of the posterior predictive distributions of the FS and RS size, respectively, as obtained by fitting the physical model described in Appendix D to the sizes shown in the figure, assuming a homogeneous (top panel) or wind-like (bottom panel) external medium, and assuming the FS to dominate at 15 GHz and the RS to dominate at 4.9 and 8.3 GHz. The shaded bands around these lines show the 68% credible interval of the posterior predictive distribution.

In the text

Fig. 7. Corner plot of the posterior probability density of the physical model parameters in the homogeneous (top panel) or wind-like (bottom panel) external medium case and assuming the FS to dominate at 15 GHz and the RS to dominate at 4.9–8.3 GHz. In each panel, the red histograms show the marginalised posterior probability densities, with black solid vertical lines showing the median and dashed lines showing the 68% credible interval or, if the latter extends to the lower (upper) extremum of the prior range, the 95% upper (lower) limit (values reported on top of the panels). The filled contours show the smallest regions containing 68% and 95% of the two-dimensional posterior probability, with the black squares showing the position of the median. The grey lines show contours of constant total jet energy assuming A = 1 cm−3. Each contour is labeled with the base-10 logarithm of the corresponding total jet energy.

In the text

Fig. B.1. Posterior probability distribution on the source size and flux density in our EVN 4.9 GHz epochs at T − T0 = 43 (RG1013 C, left), 117 (RG013 D, centre) and 261 (RG013 F, right) days. In each corner plot, the top-left and bottom-right sub-panels show histograms of the posterior samples of Fν and FWHM, with the vertical solid lines showing the median and the vertical dashed lines bracketing the 68% credible interval or, if the latter extends to 0, the 95% credible upper limit. The bottom-left sub-panel of each corner plot shows the smallest contours containing 68% and 95% of the posterior probability.

In the text

	Fig. B.2. Similar to Figure B.1, but for our EVN 8.3 GHz epochs at T − T0 = 40 days (EVN B, left) and 118 days (EVN E, right).
In the text

	Fig. B.3. Similar to Figure B.1, but for our VLBA 15 GHz epochs at T − T0 = 44 (VLBA B, upper left), 114 (VLBA C, upper right), 205 (VLBA C1, lower left) and 262 (VLBA D, lower right) days.
In the text

	Fig. C.1. Estimates of the FWHM of the GRB as a function of the total area of the synthesised beam for the EVN observations.
In the text

	Fig. D.1. Flux density (left panel) and size (right panel) of GRB 221009A and J1905+1943 in the EVN observations. The average flux density of J1905+1943 (left panel) and the 1:1 correlation (right panel) are shown as grey dashed lines for the sake of comparison.
In the text

Fig. D.2. Constraint on the viewing angle from the absence of an observed source apparent displacement in our VLBA observations. In each panel, filled contours show the displacement of the centre of the fitted Gaussian expected between 44 d and 262 d, color coded as shown in the colorbar on the right, as a function of the E/A ratio and of the off-edge viewing angle θv − θj. The red contour shows Δ(262 d)−Δ(44 d) = 0.6 mas, which represents the largest displacement compatible at 1 σ with our observations. The red contour hence contains the excluded parameter region. The upper panel row refers to a wind-like external medium, while the lower row refers to a homogeneous external medium. Each column assumes a different jet half-opening angle, given at the top of the column.

In the text