The Radio Afterglow of the Ultralong GRB 220627A

If a GRB is strongly gravitationally lensed, multiple images of the GRB would be produced (e.g., M. Oguri 2019). In the scenario where the observing instrument does not have sufficient angular resolution to resolve the lensed images, the gravitational lens signature may be superimposed onto the prompt light curve (O. M. Blaes & R. L. Webster 1992). Any flux variation from the GRB would appear at different times in the lensed images because the photons travel different path lengths. Therefore, in a spatially unresolved light curve there would be two pulses (for two lensed images), where the temporal profile of the second pulse should be similar to the first pulse, albeit weaker in intensity. The spectral properties of the two pulses in this scenario should also be similar, i.e., the lens signature is consistent across all energy bands.

The time delay between the two pulses can be explained by a combination of both the Shapiro time delay and geometric effects from differing paths traversed by the two rays; it is governed by the relation provided in L. M. Krauss & T. A. Small (1991):

where M_z = M(1 + z_L) is the redshifted lens mass (z_L is the lens redshift), r is the brightness ratio of the two lensed images, G is the gravitational constant, and c is the speed of light.

In the case of lensed GRB candidates reported to date (including GRB 220627A), the time delays between the putative lensed images are usually seconds to minutes (e.g., P. Veres et al. 2021; Y. Wang et al. 2021; X. Yang et al. 2021). In this case, using Equation (1) and assuming a range of r from 1.5 to 5 and z up to ∼9, we would infer a lens mass of ∼10⁵⁻⁷ M_⊙. We note that the optical depth of lenses in this mass range is expected to be an order of magnitude lower than those with more typical “galactic” masses (∼10¹² M_⊙) (N. Loudas et al. 2022), and that none have been conclusively identified to date. However, this suggests that lensed GRBs may provide a promising avenue to conclusively identify the first lens within this mass range (see Section 2.2 below).

Since high-energy satellites have subsecond temporal resolutions, double-pulse prompt light curves have been the basis of various claims (and disputes) of gravitationally lensed GRBs as discussed in Section 1. As an example, we show in Figure 1 the prompt light curve for GRB 220627A (the burst that is the focus of this manuscript) from Fermi/GBM. The presence of two pulses with a flux ratio r ≈ 3 (this was determined by the ratio of the pulse fluences and is typical of a doubly imaged strong lensing system; see also O. J. Roberts et al. 2022), which showed similar temporal profiles across more than one energy band, initially provided compelling evidence supporting the strong gravitational lensing scenario. However, such observations cannot rule out the possibility that the double-pulse profile may be produced by a physical emission mechanism intrinsic to the GRB (rather than by lensing), or even the possibility that the double-pulse profile does not appear in energy bands outside the observable range of the instrument (e.g., as demonstrated by GRB 220627A, the double-pulse profile was clear in the Fermi/GBM data, but was lacking in the Fermi/LAT data; Y.-Y. Huang et al. 2022).

Zoom In Zoom Out Reset image size

The smoking gun for a gravitationally lensed GRB would be the direct detection of two (or more) consistently evolving lensed images of the afterglow with the same surface brightness. The angular separation of the images is on the order of the Einstein radius:

where D_LS, D_L, and D_S are the distances between the lens and the GRB, the distance to the lens, and the distance to the GRB, respectively. For a typical event, with GRB redshift z_S = 1 and lens redshift z_L = 0.4, the Einstein radius for a 10⁶ M_⊙ point-mass lens (the typical lens mass inferred from previous claims; see above) is , where D = D_LD_LS/D_S (P. Veres et al. 2021). As a result, these systems are typically referred to as “millilens” systems. Specifically for GRB 220627A, given its redshift and that of the intervening system (which we interpret here as a potential lens), the expected Einstein radius would be θ_E ≈ 1.8 mas. Considering this typical angular scale expected from lensed GRB systems, direct detection of the two images in the prompt gamma rays or through standard X-ray/optical/radio afterglow observations would be infeasible. VLBI is currently the only way of reaching the necessary angular resolution for confirming a strongly lensed GRB event.

To confirm a strongly lensed GRB event with VLBI observations, the detection of spatially resolved components would be required. Ideally, these detections are made in at least two epochs, with at least one epoch featuring quasi-simultaneous observations at two or more frequencies. These allow for a confirmation of lensing by checking:

• 1.  
  for consistency between the spectra of two lensed images (e.g., Figure 6 of C. A. Katz et al. 1997);²¹
• 2.  
  for surface brightness preservation between the two lensed images (a necessary signature of gravitational lensing);
• 3.  
  that the two lensed images preserve their brightness ratio as they evolve and, aided by complementary low-resolution radio observations, evolve in accordance with known GRB closure relations (e.g., R. Sari et al. 1998); and
• 4.  
  that the brightness ratio of the afterglow images matches the ratio of the pulse fluences from the GRB prompt light curve (although there are certain scenarios where this may not hold, e.g., if there are significant wavelength-dependent absorption and scattering effects caused by the lens galaxy, or if the prompt emission comes from an emitting region with a large displacement from the radio emitting region—see R. Mittal et al. 2007; C. C. Cheung et al. 2014).

The combination of all these checks makes a confirmation robust against potential ambiguities that arise from using one or two checks, such as the possibility of chance spatial association with nonafterglow related sources (e.g., a background host galaxy nucleus).

The confirmation of a strongly lensed GRB candidate has scientific implications for time-delay cosmography, the physics of the lens, and our understanding of GRBs. A confirmation would provide a new example of gravitational lensing to test Einstein’s general relativity, checking for consistency in both the temporal and spatial data. Like other strong gravitationally lensed explosive transients, lensed GRBs have important applications in time-delay cosmography; the measured time delay can be used alongside a well-constrained Einstein radius to constrain the Hubble parameter H₀ (S. Refsdal 1964). This is a direct consequence of the time delay between the lensed GRB images being dependent on the mass distribution of the lens, which can be modeled from the morphology and distribution of detected components in high-resolution VLBI images, as well as geometric effects from different paths in spacetime, which is determined by the Hubble parameter—for a comprehensive review of time-delay cosmography for transients, see M. Oguri (2019).

As GRBs are found at high redshifts up to z ∼ 9 (e.g., R. Salvaterra 2015) and their prompt light curves (and any time delays) can be sampled accurately down to subsecond temporal resolution (e.g., C. Meegan et al. 2009), strongly lensed GRBs can potentially probe for lens mass and distance ranges that cannot be probed by strongly lensed quasars and supernovae (since the temporal resolution of the time delay determines the redshifted lens mass range your probe is sensitive to as shown in Equation (1)). Accompanying previous claims of lensed GRBs were suggestions of using them to study the population statistics and physics of their lenses (e.g., J. Paynter et al. 2021). In fact, a confirmation of a lensed GRB (with an observed time delay on the order of seconds to minutes) could also lead to the first confirmation of a millilens object, which implies the presence of massive compact objects, such as intermediate-mass black holes (IMBHs) and dark matter halos expected in some dark matter models (e.g., P. N. Wilkinson et al. 2001; C. Spingola et al. 2019; C. Casadio et al. 2021).

On the other hand, a conclusive rejection of a lensing explanation for a GRB with a double-pulse profile would put extra constraints on the physics of GRB central engines. For example, current GRB models that can explain multiple pulses in the prompt emission would need to be extended to explain how their emission mechanism in some GRBs can produce repeated pulses with similar temporal and spectral profiles. In the specific case of GRB 220627A, the comparison of the Fermi/GBM and Fermi/LAT data that ultimately disfavored the lensing hypothesis (Y.-Y. Huang et al. 2022) indicated that the GRB belonged to the subpopulation of ultralong GRBs.

Motivated by these broader scientific implications for confirming and characterizing a lensed GRB event, we began our radio campaign to independently test the lensing hypothesis for GRB 220627A through VLBI shortly after the GRB was detected; this campaign is described in the following section.

We observed the field of GRB 220627A with the ATCA, using the coordinates of the reported optical afterglow (S. de Wet et al. 2022) as the phase tracking center. We conducted four epochs of ATCA observations, observing the first epoch on 2022 July 5 for 10 hr (7.5 days postburst), the second epoch on 2022 July 25 for 12.5 hr (27.4 days postburst), the third epoch on 2023 September 22 for 12 hr (452.2 days postburst), and the fourth epoch on 2023 September 29 for 12 hr (459.2 days postburst), all under project C3478 (PI: J. Leung). For the first two epochs, the observations were observed at multiple bands, centered on 5.5, 9, 17, and 19 GHz, each with a bandwidth of 2 GHz, while for the final two epochs, the observations were observed only at 17 and 19 GHz, each with a bandwidth of 2 GHz. The array was in H214 configuration for the first two observations and in H168 configuration for the final two epochs.²² In these configurations, the maximum baselines are 247 and 192 m, respectively. The values do not consider baselines involving antenna CA06 (all of which exceed 4 km) since they were flagged in our data reduction process as discussed below.

We reduced the visibility data using standard Miriad procedures (R. J. Sault et al. 1995). We used PKS B1934−638 to set the flux-density scale and calibrate the bandpass response, and PKS B1313−333 (J1316−3338) to calibrate the complex gains. Prior to imaging, we flagged out all the baselines involving antenna CA06. This was because the large gap in (u, v)-coverage between these longer baselines and the other much shorter baselines in the array created difficulties in the imaging process (due to the highly irregular shape of the resulting dirty beam). We produced images using the multifrequency synthesis CLEAN algorithm (J. A. Högbom 1974; B. G. Clark 1980; R. J. Sault & M. H. Wieringa 1994) and extracted flux densities in the case of a detection by fitting a point-source model in the image plane. For the 2022 July 5 (7.5 days postburst), 2023 September 22 (452.2 days postburst), and 2023 September 29 (459.2 days postburst) observations, we did not use the 19 GHz band as it was either severely affected by poor weather conditions or radio frequency interference. However, both the 17 and 19 GHz frequency bands provided reliable data for the 2022 July 25 (27.4 days postburst) observations so we could image them together, providing a contiguous bandwidth of 4 GHz centered on approximately 18 GHz. In order to obtain a deep late-time image of the field, we also imaged the 2023 September 22 and 29 observations (mean epoch of 455.7 days postburst) at 17 GHz together, providing an effective total integration time of 15 h. The flux densities (and 3σ upper limits for nondetections) are reported in Table 1. The rms noise in all bands was higher in the first observation, compared to the second, due to poor weather conditions.

We observed the field of GRB 220627A with the MeerKAT telescope on 2022 July 6 (8.9 days postburst) and 2022 August 2 (35.7 days postburst) under project DDT-20220705-SG-01 (PI: S. Giarratana) for 2 h each. These observations were centered on 1.28 GHz (L band) with a bandwidth of 0.78 GHz. We used PKS B1934−638 for bandpass and flux-density scale calibration, and J1311−2216 for phase calibration. The archived raw visibilities were converted to the Measurement Set format with the KAT Data Access Library (katdal²³ ) for continuum imaging. We reduced the data using oxkat²⁴ (v0.3; I. Heywood 2020), where: the Common Astronomy Software Applications (CASA; J. P. McMullin et al. 2007; CASA Team et al. 2022) and TricoloUr²⁵ (B. V. Hugo et al. 2022) packages were used for measurement set splitting, cross-calibration, and flagging; CubiCal²⁶ (J. S. Kenyon et al. 2018) was used for self-calibration; and WSClean (A. R. Offringa et al. 2014) was used for continuum imaging. All processes were executed with oxkat L-band default settings. We do not detect any radio emission at the source location in both observations. The 3σ nondetection limits are given in Table 1.

We obtained observations for GRB 220627A with the VLA at three separate epochs, with 1 h of observation time per epoch, under project VLA/22B-031 (Legacy ID: AL1245; PI: J. Leung). The first epoch of observations occurred on 2022 July 13 (15.2 days postburst), at the 10 and 15 GHz (X and Ku) frequency bands, with the array in A→D configuration.²⁷ We then obtained a second epoch at the 6 GHz (C band) frequency band on 2022 August 5 (39.0 days postburst) and a third epoch at the 15 GHz (Ku band) frequency band on 2022 August 11 (45.1 days postburst); the array was in D configuration for both the second and third epoch of observations. For each observation, we used 3C 286 to calibrate the bandpass response and flux-density scale; we also used PKS B1313−333 (J1316−3338) to calibrate the complex gains.

Prior to imaging, we obtained the pipeline-calibrated²⁸ visibility data from the data archives (see Data Availability section). We used CASA to verify the quality of the pipeline-calibrated data and also to image the field using the wide-field, multifrequency synthesis CLEAN algorithm. For the 2022 July 13 (15.2 days postburst) observation, due to problems associated with the unbalanced (u, v) coverage of the hybrid array, we were only able to obtain reliable data by using just the antennas that had already moved to the D configuration (i.e., we flagged all antennas located more than 2 km from the center of the array). After further accounting for antennas under maintenance and antennas flagged by the automated pipeline, the imaging process for the 10 and 15 GHz observations on 2022 July 13 (15.2 days postburst) used data from only 12 and 13 antennas, respectively; for the 2022 August 5 (39.0 days postburst) observation (at 6 GHz), the imaging process used data from 21 antennas; and for the 2022 August 11 observation (45.1 days postburst) at 15 GHz, the imaging process used data from 22 antennas. Further self-calibration was not performed for all observations since weak field sources in addition to the presence of poorly modeled extended sources led to limited improvements in the gain solutions. Due to the presence of imaging artifacts caused by the sidelobes of a bright, extended source outside the primary beam, PKS J1325−3236, we phase rotated the 2022 August 5 (6 GHz band) visibility data to J2000 13:25:32.8, −32:30:45 before imaging so we could deconvolve both the target and PKS J1325−3236 in our CLEAN cycles. For each observation, we imaged using the full bandwidth, i.e., 4 GHz for the observations centered on 6 and 10 GHz and 6 GHz for the observations centered on 15 GHz. In the case of a detection, we extracted the flux densities by fitting a point-source model in the image plane and report them in Table 1 (for nondetections, we give the 3σ upper limit).

We obtained high-resolution observations for GRB 220627A with the VLBA on two consecutive days under project VLBA/22B-032 (Legacy ID: BL296; PI: J. Leung). These observations were 6 hr in duration, starting at 22:30 UT on both 2022 July 8 and 9 (corresponding to a mean epoch of 11.7 days postburst). All standard VLBA antennas,²⁹ with the exception of the Kitt Peak station, which was out due to wildfires, were available for both observations. The longest baseline was from Maunakea to Saint Croix, giving a projected baseline of ∼8600 km or ∼430 Mλ. For these observations, the VLBA was configured to use the digital downconverter (DDC) observing system, observing at the central frequency of 15.168 GHz (Ku band), with a total bandwidth of 512 MHz (split into four 128 MHz subbands) in two circular polarizations (right- and left-circular polarizations). Using a 2 bit sampler, the resulting data rate was 4 Gbps. The data were correlated after the observations using DiFX (A. T. Deller et al. 2007) with an integration time of 1 s and a spectral resolution of 1 MHz at Socorro, New Mexico, USA.

We carried out standard reduction procedures on the correlated datasets from both observations simultaneously using AIPS (E. W. Greisen 2003). We applied standard a priori amplitude calibration using the provided system temperature and gain curve tables. The instrumental delays were determined using the fringe finders J0555+3948 and J1642+3948. The latter was also used for solving the bandpass. Since the expected flux density of our target, GRB 220627A, was <0.5 mJy as informed by our earlier ATCA observations, we phase-referenced our experiment using the calibrator J1324−3235, which was observed in 45 s scans as part of 100 s duty cycle. Although this calibrator had a low unresolved flux-density measurement at 8.6 GHz (50 mJy), it was chosen due to its close separation (032) from the target position, which was particularly important given the low elevation of the source with respect to the VLBA antenna stations. To attain a better signal-to-noise ratio, we combined all the subbands during our global fringe fitting process on the fringe finders and phase-reference calibrators in order to solve for the residual phases, delays, and rates. The solutions from the phase-reference calibrator were interpolated to both the target and our check source, PKS B1313−333 (J1316−3338), located 175 away from the target. Finally, we used Difmap (M. C. Shepherd 1997) with a natural weighting scheme (after trying several other robust weighting schemes to search for sources in the field of interest) to obtain the final images of both the target and the check source. The resulting image had a beam size of 0.49 mas × 1.51 mas with a position angle of −54.

In our VLBA observations, we were unable to identify any sources that were related to the GRB 220627A event within the optical³⁰ and radio³¹ positional uncertainties. The final target image achieved an rms sensitivity of 51 μJy beam⁻¹ and no noise peaks exceeded the 5σ level (i.e., ∼250 μJy beam⁻¹). After performing self-calibration iterations on the check source (located further from the phase-reference calibrator than the target), we found that the phase errors were no more than ∼10° for most antennas, suggesting that possible problems with phase transfer were not responsible for the nondetection. The 5σ upper limit is reported in Table 1. We verified this result in an independent calibration of the data (as a blind test) using standard reduction procedures in CASA; we attained a similar rms sensitivity (and upper limit), arriving at the same conclusions.

GRB 220627A was originally proposed as a candidate lensed GRB based on similar pulse profiles observed in the Fermi/GBM data (O. J. Roberts et al. 2022). While this interpretation was later disfavored after assessing the spectra inferred for both pulses using data from Fermi/LAT (Y.-Y. Huang et al. 2022), we check to see whether a similar conclusion could have been drawn from the Fermi/GBM data alone.

In particular, we analyze the Fermi/GBM data of both the NaI#3 and BGO#0 detectors. We used the Fermi Science Tools v.2.2.0 to analyze the Fermi/GBM data following the standard procedure³² (e.g., M. E. Ravasio et al. 2019, 2024). We extracted the light curves of NaI#3 in the 10–900 keV energy range and of BG#0 in the 300–30,000 keV energy range. The background-subtracted light curve of these two detectors is shown in the top and bottom panels, respectively, of Figure 1. The prompt emission shows two distinct emission episodes, each lasting a few hundred seconds, separated by a much longer period of quiescence. The two emission episodes have quite different temporal structures, with the first one being more variable and the second one apparently composed of two subpulses.

We extracted the time integrated spectrum of the two emission episodes by selecting two time intervals, [−6, 300] and [910, 1096] s, shown by the shaded regions in Figure 1. Given the large separation of the two time intervals, the background spectrum was estimated by considering four intervals bracketing the two emission episodes: [−450, −200], [400, 500], [700, 850], and [1200, 1300] s. A third-order polynomial was fit to the spectra contained in these time intervals and then interpolated to the time region selected for the spectral analysis.

We performed the spectral analysis with Xspec v.12.14.0 (K. A. Arnaud 1996). We considered the energy range between 10 and 900 keV for the NaI spectrum and excluded the 30–40 keV range due to the contamination by the Iodine K-edge at such energies. The energy range 300 keV–30 MeV was considered for the BGO spectrum. We allowed for a free intercalibration constant factor between the NaI and BGO spectra. We adopted the PG-stat, i.e., Poisson statistic for the source counts and Gaussian statistic for the background counts, in performing the fit.

The spectral model adopted is a smoothly broken power law:

in photons cm⁻² s⁻¹ keV⁻¹ units. α and β are the spectral index of the power laws below and above, respectively, the characteristic break energy E_break. The model parameter n sets the curvature between the two power laws. We assume n = 2, which resembles the curvature of the typical Band function, which has been fitted to GRB spectra (D. Band et al. 1993, see also M. E. Ravasio et al. 2018).

The normalization constant in Equation (3) is defined as:

where L is the luminosity in units of erg s⁻¹ in the 1 keV−10 MeV rest-frame energy range and d_L(z) is the luminosity distance corresponding to the redshift z of the source. Therefore, the luminosity L is one of the four free parameters of the spectral model.

The peak energy of the spectral energy distribution (SED), i.e., the energy where E² × dN/dE peaks, is:

We report the fit results for the two time intervals in Table 2. The posterior distributions of the fit parameters are shown in Figure 2 and the resulting fits are shown in Figure 3. In addition to the aforementioned differences in temporal structures between the first and second episodes (as shown in Figure 1), the spectra taken at the two time intervals are clearly distinct from each other, e.g., the α slope parameter differs between the two intervals with a significance of >5σ. It is likely that the initial lensing hypothesis proposed in O. J. Roberts et al. (2022) was based off the analysis of the Fermi/GBM data with lower-resolution binning, where the temporal structures between the two episodes are more similar. Therefore, based on the reanalysis of the Fermi/GBM data alone—without requiring the additional information from Fermi/LAT as in Y.-Y. Huang et al. (2022)—we can reject the lensed GRB hypothesis and verify the ultralong GRB classification.

Table 2. Spectral Parameters of the Fit of the Prompt Emission Spectrum of the Two Emission Episodes

Interval	log(L)	Epeak	α	β	factor	PG-stat(dof)
(s)	(erg s−1)	(keV)
[−6, 300]						704(219)
[910, 1096]						337(219)

Note. Errors represent the 95% confidence interval for each parameter value; 'dof' is degrees of freedom.

Download table as:  ASCIITypeset image

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In addition to informing our VLBI observations, our broadband radio light curves provide a detailed view of the afterglow of GRB 220627A between ∼4 and 450 days postburst. Here we assess implications of these observations on the properties of the afterglow and the progenitors of ultralong GRBs.

4.2.1. Basic Properties of the Radio Afterglow

Radio frequency light curves and spectra of the GRB 220627A afterglow at 7.5, 15, and 27 days postburst are shown in Figure 4. Our ATCA observations 7.5 post-trigger (light blue spectrum in the right panel of Figure 4) indicate a positive spectral slope α > 1.5 between 9 and 17 GHz, where F_ν ∝ ν^α and we quote a 3σ lower limit. Assuming synchrotron emission, this is a strong indication of self-absorption at these frequencies and at this epoch, with the self-absorption frequency ν_a > 17 GHz. At later times, our observations still indicate a rising spectrum, but the available limits can also accommodate a shallower slope of F_ν ∝ ν^1/3 (Figure 4, right). Due to a lack of detections, we cannot constrain the temporal evolution at ν < 15 GHz. The 15–17 GHz light curves, on the other hand, reveal an unexpected flattening in the decay of the radio flux density at t ≳ 15 days (Figure 4, left): while the evolution between 7.5 and 15 days is consistent with a steep F_ν ∝ t⁻² decay, the subsequent detections indicate a shallower F_ν ∝ t^−1/2.

Zoom In Zoom Out Reset image size

4.2.2. Comparison to the Afterglow Model of S. de Wet et al. (2023)

Previously, S. de Wet et al. (2023, DW23 hereafter) presented a best-fit model for the afterglow of GRB 220627A. Specifically, they used the ScaleFit software package (G. Ryan et al. 2015) to model a set of X-ray, optical, and radio afterglow data (their radio data points are a preliminary subset of what we present in this paper). Overall, their modeling favored a homogeneous circumburst medium and converged to the following afterglow model parameters: slope of the electron energy distribution , the isotropic-equivalent kinetic energy of the outflow erg, circumburst medium density cm⁻³, the microphysical shock parameters relating to the fractional internal energy in accelerated electrons and magnetic fields and ×10⁻³, jet half-opening angle , and host galaxy extinction A_V,host = 0.09 ± 0.03 mag.

In Figure 5, we plot the radio data from Table 1, optical data from DW23, and the X-ray data from the Swift Burst Analyzer³³ (P. A. Evans et al. 2010). Also shown are model light curves and spectra obtained with BOXFIT³⁴ v. 2.0 (H. van Eerten et al. 2012), with the same parameters as the best-fit model of DW23.

Zoom In Zoom Out Reset image size

Our 7.5 days radio spectrum (epoch 4 in Figure 5) is in tension with the predictions from the DW23 model, which overpredicts the flux density at 10 GHz by at least a factor of ∼2 at this time. Moreover, the temporal evolution revealed by our 15–17 GHz observations does not match that predicted by the model: the relatively early jet break in the DW23 model causes the radio light curve to drop steeply as F_ν ∝ t^−p (where p = 2.0 in their model) after ∼5–7 days (corresponding to the time in their model after which both the injection frequency, ν_m, and the self-absorption frequency, ν_a, cross the observing band). While this is consistent with the observed decay between 7.5 and 15 days, the subsequent flatter evolution remains unexplained. We note that we found a similar disagreement also with the alternative model presented in Appendix A of DW23, which assumes a wind-like external medium.

Given that the flux density at 17 GHz is seen to vary, with a decay at least as steep as F_ν ∝ t^−1/2 from 27 to 445 days (see Figure 4, left), and given the positive spectral slope implied by our upper limits at lower frequency, we can exclude that the excess emission with respect to the DW23 model is due to star formation in the host galaxy.

In addition, at least two of our upper limits below 17 GHz (at 7.5 and 15 days) are inconsistent with the DW23 model. However, given the good match between the model and higher-frequency observations at these epochs, a possible explanation of such an inconsistency could be traced back to the simplified description of the electron population in the afterglow model: the energy distribution of electrons in GRB afterglow models (e.g., R. Sari et al. 1998; A. Panaitescu & P. Kumar 2000; H. van Eerten et al. 2012) is described as a nonthermal power law, but physically we expect a lower-energy thermal population to be present in these relativistic outflows (e.g., D. Giannios & A. Spitkovsky 2009; L. Sironi et al. 2013; D. C. Warren et al. 2017; H.-X. Gao et al. 2024; L. Rhodes et al. 2025). Even in the case in which the latter electrons constitute a nondominant fraction over the total, with a negligible contribution to the radiated power, they could still affect the low-energy spectrum by increasing the self-absorption optical depth (D. Eichler & E. Waxman 2005; S. M. Ressler & T. Laskar 2017; D. C. Warren et al. 2018, 2022; B. Margalit & E. Quataert 2021, 2024). For that reason, we do not consider these inconsistencies as severe.

4.2.3. Explaining the Shallow Radio Decay

We discuss below three scenarios that can explain the shallow decay observed in radio at t > 15 days. For each scenario, we describe the modeling choices and the implications on our understanding of the source.

Energy injection from an additional, slower ejecta component. The decay of the radio light curve could become shallower due to the energy injection occurring when a slower ejecta component catches up to the external shock (A. Panaitescu et al. 1998; R. Sari & P. Mészáros 2000). First, we can estimate the Lorentz factor of such ejecta for the collision to happen at an observed time t_coll ∼ 15 days, assuming a similar ejection time as the main ejecta that caused the peak of the afterglow emission. Assuming the dynamics of the blast wave to be described, before the collision, by the R. D. Blandford & C. F. McKee (1976) self-similar solution in a homogeneous external medium, at an observer time t the part of the blast wave closest to the observer is located (R. Sari et al. 1998) at a radius from the central engine, where m_p here is the proton mass. The radius of the slower ejecta component is instead . Equating these two radii and setting t = t_coll, we obtain the required ejecta Lorentz factor:

In words, ejecta as fast as Γ_sl−ej ∼ 10 catch up with the blast wave at as late as t ∼ 15 days. We describe the energy injection as an increase in the isotropic-equivalent blast wave energy with observer time³⁵ at t > t_coll, parameterized as . This modifies the blast wave Lorentz factor evolution such that the shock downstream Lorentz factor scales as Γ(t) ∝ t^(4δ−3)/8 and the blast wave radius as r ∝ t^(δ+1)/4. In the post-jet-break phase, neglecting sideways expansion of the shock and focusing on the ν_a, ν_m < ν < ν_c spectral regime, the flux density is proportional to F_ν ∝ Γ^2(p+1)r³, which leads to a light-curve scaling F_ν ∝ t^{((7+4p)δ−3p)/4}. With p ∼ 2, a scaling F_ν ∝ t^−1/2 can be reproduced by setting δ ∼ 4/15 ≈ 0.27. This scaling is expected (A. Panaitescu et al. 1998; R. Sari & P. Mészáros 2000; B. Zhang et al. 2006) if the slow ejecta features a velocity profile such that its mass above a Lorentz factor Γ is M( > Γ) ∝ Γ^−s with s = (7δ + 3)/(3–δ) ∼ 73/41 ≈ 1.8. Since the observer time increases fourfold during the shallow decay phase, the energy increase due to injection is just ∼4^4/15 ≈ 145%.

These calculations show that a relatively slow shell ejected shortly after the main ejecta and containing about half as much energy could produce the observed flattening in the radio light curve. The presence of a tail of slower ejecta carrying a substantial energy is a common expectation from the unsteady central engines behind leading GRB jet prompt emission scenarios, such as the internal shocks scenario (M. J. Rees & P. Meszaros 1994) and the internal-collision-induced magnetic reconnection and turbulence (ICMART; B. Zhang & H. Yan 2011) scenario. Unfortunately, model predictions are not yet sufficiently detailed as to answer the question whether the requirements found here are easily satisfied in a typical central engine.

Changing external density profile. The flattening in the radio light-curve decay could stem from a steepening in the external density profile, which could signal, for example, the fact that the blast wave is emerging from the molecular cloud where the progenitor has formed. Assuming the external density to scale as n ∝ r^−k, the Lorentz factor in the shock downstream then evolves (R. D. Blandford & C. F. McKee 1976) as Γ ∝ r^−(3−k)/2. In terms of observer-frame time, we then have r ∝ t^1/(4−k) and Γ ∝ t^{−(3−k)/2(4+k)}. In the post-jet-break phase the flux density evolves as F_ν ∝ Γ^2(p+1)r^3−k (the 3 − k exponent here reflects the slower increase in the number of shocked electrons with radius when k > 0), so that F_ν ∝ t^{−(3−k)p/(4−k)}. This shows that the t^−1/2 decay can be reproduced, assuming p ∼ 2, if the external density decreases with a slope k = (6p − 4)/(2p − 1) ∼ 8/3.

An additional, wider ejecta component. As an alternative scenario, we consider the possibility of a wider and slower component in the ejecta, on top of the fast and narrow component whose afterglow emission is described by the DW23 model. These ejecta decelerate at a larger radius with respect to the narrow jet, enhancing the observed radio flux density. We model the dynamics and emission of such an additional component following A. Panaitescu & P. Kumar (2000), including the ejecta coasting phase, the transrelativistic deceleration of the blast wave, the effect of synchrotron self-absorption on the emission (but still neglecting thermal electrons), and photon propagation effects. In Figure 6, we show the model light curve obtained by adding the flux densities predicted by such a model to those of the DW23 model, keeping the microphysical parameters fixed to values similar³⁶ to DW23, that is, _e = 0.076, _B = 4.6 × 10⁻³, and p = 2.02. The wide jet component is assigned an isotropic-equivalent energy E_iso,w = 3 × 10⁵³ erg, an initial bulk Lorentz factor Γ_0,w = 3, and a half-opening angle θ_j,w = 15°. Its collimation-corrected energy is therefore erg. For comparison, the narrow DW23 jet has E_iso = 9 × 10⁵³ erg and θ_j = 4.5° (and a formally infinite initial bulk Lorentz factor), hence a collimation-corrected energy of erg. The addition of the wide component therefore increases the energy budget by about a factor of 4. The inclusion of this wider component leads to a clear improvement in the modeling of the late flattening in the radio light curves, as shown in the comparisons of the residuals in Figures 5 and 6.

Zoom In Zoom Out Reset image size

The presence of a slow, wide ejecta component that carries a similar (or even somewhat larger) collimation-corrected energy as the inner, narrow jet is expected in the case where the jet spends a large fraction of its total energy in breaking out of the progenitor star (C. D. Matzner 2003; O. Bromberg et al. 2011; R. Harrison et al. 2018; O. Gottlieb et al. 2021). In such a scenario, the energy is channeled to a hot cocoon that shrouds the jet during the propagation. After the jet breakout, the cocoon blows out producing a slow, wide ejecta component (O. S. Salafia & G. Ghirlanda 2022), which could explain the scenario behind our proposed modeling. For the cocoon to contain a large fraction of the total energy, the jet head must have propagated at a Newtonian velocity for most of its journey through the progenitor envelope (C. D. Matzner 2003; O. Bromberg et al. 2011; O. Gottlieb et al. 2021; O. S. Salafia & G. Ghirlanda 2022). In such a regime, the time for the jet head to break out can be estimated based on the numerically calibrated model of R. Harrison et al. (2018, their Equation (16)), which gives a breakout time (relative to the birth time of the jet):

where E_tot,52 is the total energy in the ejecta (narrow plus wide) in units of 10⁵² erg, T₃ is the total jet duration in units of 1000 s, is the half-opening angle of the jet at its base (i.e., at launch, before propagating through the stellar envelope) in units of 10°, R₁₁ is the stellar radius in units of 10¹¹ cm, and M₁₀ is the stellar mass in units of 10 M_⊙.

In this regime, the energy in the cocoon (and hence in the wide ejecta) is approximately E_w ∼ E_tott_bkt/T. The total jet duration T can be estimated as the sum of the breakout time plus the GRB duration, T ∼ t_bkt + T₉₀/(1 + z). Using E_tot = E + E_w, these relations yield:

and

Here we used E_w and E from the proposed modeling and T₉₀ ∼ 1100 s from DW23, and we kept the initial half-opening angle at 10° (the jet is then expected to be collimated during propagation; see, e.g., O. Bromberg et al. 2011, O. S. Salafia & G. Ghirlanda 2022). If indeed the wide component in our proposed modeling originates from a cocoon, the long jet duration and large stellar radius obtained from these arguments are then consistent with expectations for the blue supergiant progenitor scenario proposed for ultralong GRBs (e.g., B. Gendre et al. 2013; K. Ioka et al. 2016; R. Perna et al. 2018). In particular, we note that typical Wolf–Rayet stars have radii of only ∼1 R_⊙ compared to tens to a few hundreds of solar radii for blue supergiants. This is similar to the conclusion reached by L. Piro et al. (2014) in their study of ultralong GRB 130925A (with gamma-ray activity lasting for at least ∼7 ks); they found evidence for the presence of a hot cocoon that supported the blue supergiant progenitor scenario. More recently, by conducting simulations of evolutionary pathways for different possible progenitor scenarios, the findings in A. K. Ror et al. (2024) similarly support blue supergiant stars as a plausible progenitor for ultralong GRBs, including GRB 221009A, which was the focus of their work.

The aim of our VLBA experiment was to provide a test of the lensing scenario for GRB 220627A that was independent of the conclusions drawn from the analysis of the high-energy spectra (see Section 4.1; see also Y.-Y. Huang et al. 2022). From the light-curve evolution shown in Figure 6 (left), we estimate the unresolved flux density at the time and frequency of our VLBA observation (15 GHz, 11.7 days postburst) to be ∼150 μJy. This suggests that in the ultralong GRB (nonlensing) scenario, we should expect a compact radio source at ∼150 μJy, and in the lensing scenario two compact sources at even dimmer flux-density levels of ∼40 and ∼110 μJy (as determined by their gamma-ray fluence brightness ratio r ≈ 3).

Our VLBA observations yielded a nondetection down to a 5σ limit of ∼0.3 mJy beam.⁻¹. Within the uncertainty region of lower-resolution optical and radio detections (see Section 3.4), there were a few noise peaks exceeding the 3σ level, but none above the 5σ level. This means we were unable to claim a detection or a nondetection of an afterglow source (or any source) at the ∼0.3 mJy level; the lensing and ultralong GRB scenarios could not be distinguished, so our VLBI experiment yielded a nonconclusive result.

This experiment was the first time VLBI was used to search for the multiple images of a candidate strongly gravitationally lensed GRB afterglow. The only lensed GRB candidate that had VLBI follow-up prior to this event was GRB 200716C (S. Giarratana et al. 2023). For that event, as the afterglow had faded by the time of the observation, the European VLBI Network (EVN) and e-Merlin were used to search the candidate host galaxy for a possible IMBH that could act as the lens, rather than directly for the multiple images of the afterglow. Their nondetection (5σ upper limit of <50 μJy beam⁻¹) also yielded an inconclusive result, unable to confirm the presence of a possible IMBH lens in the candidate host galaxy. These results highlight the difficulty of such an experiment and the need for improvements in future programs. Below, we evaluate the outcomes of our VLBI campaign for GRB 220627A and discuss the possible improvements that can be implemented in the areas of angular resolution, sensitivity, and trigger timing (including real-time correlation capabilities).

4.3.1. Angular Resolution: Probing Different Lens Configurations

Our resulting VLBA image had a synthesized beam size of 0.49 mas × 1.51 mas. In Figure 7, we show the expected Einstein radius for GRB 220627A in the lensing scenario as a function of the lens redshift, calculated using Equations (1) and (2), compared to our major-axis beam size. With the major-axis beam size being smaller than the expected Einstein radius out to z ≈ 2.8 (n.b., z_GRB = 3.084), our VLBA observations had sufficient angular resolution to resolve the lensed images for much of the lens redshift parameter space. As an intervening absorber was identified at z = 2.665 via the presence of S iv, S ii, C iv, Fe ii, and Al ii lines in the optical spectrum (L. Izzo et al. 2022; S. de Wet et al. 2023), we considered this system the likely location of the lens during our radio campaign. At this redshift, the Einstein radius would be θ_E ≈ 1.8 mas, so the two distinct lensed images of the afterglow (in the lensing scenario) would have been resolved by our VLBA observations. However, in the extreme scenarios, where the lens redshift is greater than z_L ≈ 2.8, the lensed images may be unresolved. In this scenario—if detected—the lensed event could potentially still be confirmed by modeling a multicomponent source, but this may be at the expense of accuracy in both the astrometry and source structure, which are critical for modeling of the lens (and by extension, time-delay cosmography).

Zoom In Zoom Out Reset image size

In Figure 7 we also plot 100 possible Einstein radii functions (i.e., as Einstein radius as a function of lens redshift) for a typical GRB at redshift z = 1.5 (the redshift where star formation and long GRB production is expected to be dominant). Each different Einstein radius function is produced by assuming different, random combinations of brightness ratio (selected from 1.1 < r < 5.0) and time delay (selected from 0.6 s < Δt < 1000 s) for the z = 1.5 GRB; the range chosen for these two variables reflects the observed range of values in previous claims of lensed GRB candidates. It is clear from this plot that the specific angular resolution required to test the lensing hypothesis for future candidates will be quite different depending on the lens configuration; nonetheless, for a large portion of such configurations, the milliarcsecond scales probed by VLBI will be very useful for testing the lensing hypothesis.

Future improvements to the program should include the optionality to add longer baselines to the array, thereby improving the resolution; for example, adding EVN telescopes like Effelsberg to the VLBA to improve east–west (u, v) coverage, or similarly, adding the East Asia VLBI Array to the Long Baseline Array (LBA) in Australia to improve north–south (u, v) coverage.

4.3.2. Sensitivity: Improvements with Current and Future Arrays

4.3.3. Trigger Timing: The Need for Sensitive Real-time Capabilities