Heavy-element production in a compact object merger observed by JWST

The mergers of binary compact objects such as neutron stars and black holes are of central interest to several areas of astrophysics, including as the progenitors of gamma-ray bursts (GRBs)1, sources of high-frequency gravitational waves (GWs)2 and likely production sites for heavy-element nucleosynthesis by means of rapid neutron capture (the r-process)3. Here we present observations of the exceptionally bright GRB 230307A. We show that GRB 230307A belongs to the class of long-duration GRBs associated with compact object mergers4–6 and contains a kilonova similar to AT2017gfo, associated with the GW merger GW170817 (refs. 7–12). We obtained James Webb Space Telescope (JWST) mid-infrared imaging and spectroscopy 29 and 61 days after the burst. The spectroscopy shows an emission line at 2.15 microns, which we interpret as tellurium (atomic mass A = 130) and a very red source, emitting most of its light in the mid-infrared owing to the production of lanthanides. These observations demonstrate that nucleosynthesis in GRBs can create r-process elements across a broad atomic mass range and play a central role in heavy-element nucleosynthesis across the Universe.

A similar emission-like feature is also visible in the later epochs of X-shooter observations of AT2017gfo [23], measured at 2.1 microns by [37]. This further strengthens both the kilonova interpretation and the redshift measurement of GRB 230307A ( Figure 3). We interpret this feature as arising from the forbidden [Te III] transition between the ground level and the first fine structure level of tellurium, with an experimentally-determined wavelength of 2.1044 microns [38]. The presence of tellurium is plausible, as it lies at the second peak in the r-process abundance pattern, which occurs at atomic masses around A ≈ 130 [39]. It should therefore be abundantly produced in kilonovae, as is seen in hydrodynamical simulations of binary neutron star mergers with nucleosynthetic compositions similar to those favoured for AT2017gfo [40]. Furthermore, the typical ionization state of Te in kilonova ejecta is expected to be Te III at this epoch because of the efficient radioactive ionization [41]. Tellurium has recently been suggested as the origin of the same feature in the spectrum of AT2017gfo [42]. A previous study [37] also identified this tellurium transition and noted that the observed feature is most likely two blended emission lines. However, alternative transitions from heavy r-process elements have been considered for this feature [e.g. Ce III ; 37]. Tellurium can also be produced via the slower capture of neutrons in the sprocess. Indeed, this line is also seen in planetary nebulae [43]. The detection of [Te III] 2.1 µm extends on the earlier detection of strontium, a first r-process peak element, in the early time photospheric emission of AT2017gfo [44]. The mass of Te III estimated from the observed line flux is ∼ 10 −3 M (see Supplementary Information 6.2). Although weaker, we also note that the spectral feature visible at 4.5 microns is approximately consistent with the expected location of the first peak element selenium and the near-third peak element tungsten [45].
Detailed spectral fitting at late epochs is challenging because of the breakdown of the assumptions regarding local thermodynamic equilibrium (LTE) which are used to predict kilonova spectra at earlier ages, as well as fundamental uncertainties in the atomic physics and electron transitions in the highly complex electron shells of r-process elements. However, these observations provide a calibration sample for informing future models. The red continuum emission indicates a large opacity in the mid-IR at low temperatures, e.g., ∼ 10 cm 2 g −1 at ∼ 700 K. Since the Planck mean opacity under LTE is expected to be less than 1 cm 2 g −1 at 1000 K [46], this large opacity may suggest that lanthanides are abundant in the ejecta. Indeed, it has been shown that non-LTE effects can increase the lanthanide opacity in mid-IR at late times [47]. Therefore, systematic studies of non-LTE opacity are necessary to answer the question whether lanthanides are the origin of the red emission at this epoch.
A fit to our combined MUSE + JWST data for the host galaxy (Supplementary Information) suggests a relatively low mass (∼ 2.5 × 10 9 M ) dominated by an older ∼ 10 Gyr population, but with a second more recent burst of star formation. These properties are entirely consistent with the population of short GRB host galaxies [48,49]. The host normalised offset of the burst from the host galaxy places it in the top 10% of those seen for short GRBs [48,50]. The offset could readily be achieved by a binary with a velocity of a few hundred km s −1 and a merger time of hundreds of millions of years. Alternatively, in the second epoch of JWST observations, a faint source is detected in the F150W images at the location of the transient. This may represent continued emission from the transient. However, its absolute magnitude of M F 150W ∼ −8.5 is comparable to the absolute magnitude of globular clusters in which dynamical interactions could be at play to create merging systems at enhanced rates [51]. Future observations should readily distinguish between a fading afterglow or underlying cluster.
It is striking that GRB 230307A is an extremely bright GRB, with only the exceptional GRB 221009A being brighter [52]. Of the ten most fluent Fermi/GBM GRBs, two are now associated with kilonovae (230307A and 211211A), three are associated with supernovae, and the nature of the remaining five appears unclear (see Table 4). For bright GRBs, there may be a significant contribution from mergers. Indeed, such a conclusion can also be reached by considering the energetics. Both GRB 230307A and GRB 211211A have isotropic equivalent energies of E iso > 10 51 erg. The majority of local GRBs for which the connection between GRBs and core-collapse supernovae has been established are much less energetic (typically E iso > 10 49−50 erg) and it has been suggested that they represent a separate population powered via different processes [53]. For more energetic bursts in the local Universe (where supernovae can still readily be discovered) the fraction of long GRBs with and without supernovae appear similar (see Supplementary Information). If a substantial number of long GRBs are associated with compact object mergers, they provide an essential complement to gravitationalwave (GW) detections. Firstly, joint GW-GRB detections, including long GRBs, can push the effective horizons of GW detectors to greater distances and provide much smaller localisations. In the case of GRB 230307A, the distance of 300 Mpc could have provided a robust detection in gravitational waves for the relevant O4 sensitivity [17,54]. Secondly, long GRBs can be detected without GW detectors, providing a valuable route to enhancing the number of kilonova detections. Thirdly, JWST can detect kilonova emission at redshifts substantially beyond the horizons of the current generation of GW detectors, enabling the study of kilonovae across a greater volume of the Universe.
The duration of the prompt γ-ray emission in these mergers remains a challenge to explain. In particular, the natural timescales for emission in compact object mergers are much shorter than the measured duration of GRB 230307A. Previously suggested models that may also explain GRB 230307A include magnetars [55,56], black hole -neutron star mergers [57,58], or even neutron star -white dwarf systems [19,59]. Recent results have also shown that the jet timescale does not directly track the accretion timescale in compact object mergers, and that long GRBs may be created from very short lived engines [60], and hence from binary neutron star mergers without magnetars.
There is evidence that the kilonova in GRB 230307A produced elements across a wide range of atomic mass. The detection of second peak elements in the spectrum of a kilonova demonstrates that nuclei with atomic masses around A ∼ 130 are being created in the mergers of compact objects. Many second peak elements have important biological roles. For example iodine is essential for mammals and may have been used by the single cell Last Universal Common Ancestor [61]. The creation of these elements in compact object mergers, which can have long delay times, may have important consequences for the time at which certain evolutionary channels become plausible. The burst begins very hard, with the count rate dominated by photons in the hardest (100-900 keV) band, but rapidly softens, with the count rate in the hard band being progressively overtaken by softer bands (e.g. 8-25 keV and 25-100 keV) beyond ∼ 20s. This strong hard-to-soft evolution is reminiscent of GRB 211211A [20] and is caused by the motion of two spectral breaks through the γ-ray regime (see Methods). The right panel shows the X-ray light curves of GRBs from the Swift X-ray telescope. These have been divided by the prompt fluence of the burst, which broadly scales with the X-ray light curve luminosity, resulting in a modest spread of afterglows. The greyscale background represents the ensemble of long GRBs. GRB 230307A is an extreme outlier of the > 1000 Swift-GRBs, with an extremely faint afterglow for the brightness of its prompt emission. Other merger GRBs from long bursts occupy a similar region of parameter space. This suggests the prompt to afterglow fluence could be a valuable tool in distinguishing long GRBs from mergers and those from supernovae. The putative host is the bright face-on spiral galaxy, while the afterglow appears at a 30-arcsecond offset, within the white box. The lower panel shows cut-outs of the NIRCAM data around the GRB afterglow location. The source is faint and barely detected in the bluer bands but very bright and well detected in the red. In the red bands, a faint galaxy is present northeast of the transient position. This galaxy has a redshift of z = 3.87, but we consider it to be a background object unrelated to the GRB (see Supplementary Information).   Supplementary Information). The decay rates in both the optical and IR are very similar to those in AT2017gfo. These are too rapid for any plausible afterglow model (e.g. as a power-law, they decay faster than t −3.5 over a prolonged period). There is also good agreement in the late time slope between the measurements made at 4.4 microns with JWST and 4.5 microns for AT2017gfo with Spitzer [31]. 13 1 Methods

Observations
Below we outline the observational data that were used in this paper. Magnitudes are given in the AB system unless stated otherwise. We utilize cosmology resulting from the Planck observations [62]. All uncertainties are given at the 1σ level unless explicitly stated.

Fermi/GBM data analysis
In Figure 1, we plot the light curve of GRB 230307A as seen by the Fermi/GBM in several bands, built by selecting Time Tagged Event (TTE) data, binned with a time resolution of 64 ms. The highlighted time interval of 3-7 s after trigger are affected by data loss due to the bandwidth limit for TTE data [73].
For the spectral analysis, we made use of the CSPEC data, which have 1024 ms time resolution. Data files were obtained from the online archive 1 . Following the suggestion reported by the Fermi Collaboration [73], we analysed the data detected by NaI 10 and BGO 1, which had a source viewing angle less than 60 • , and excluded the time intervals affected by pulse pile-up issues (from 2.5 s to 7.5 s). The data extraction was performed with the public software gtburst, while data were analysed with Xspec. The background, whose time intervals have been selected before and after the source, was modelled with a polynomial function whose order is automatically found by gtburst and manually checked. In the fitting procedure, we used intercalibration factors among the detectors, scaled to the only NaI analysed and free to vary within 30%. We used the PG-Statistic, valid for Poisson data with a Gaussian background. The best-fit parameters and their uncertainties were estimated through a Markov Chain Monte Carlo (MCMC) approach. We selected the time intervals before and after the excluded period of 2.5-7.5 s due to instrumental effects. In particular, we extracted 2 time intervals from 0 to 2.5 s (1.25 s each) and 14 time intervals from 7.5 s to 40.5 (bin width of 2 s, except the last two with integration of 5 s to increase the signal-to-noise ratio), for a total of 16 time intervals. We fitted the corresponding spectra with the two smoothly broken power-law (2SBPL) function [74,75], which has been shown to successfully model the synchrotron-like spectral shape of bright long GRBs, including the merger-driven GRB 211211A [20].
From our spectral analysis we found that all spectra up to ∼ 20 s are well modelled by the 2SBPL function, namely they are described by the presence of two spectral breaks inside the GBM band . In particular, in the time intervals between 7.5 and 19.5 s, the low-energy break E break is coherently decreasing from 304.3 +5.2 −2.6 keV down to 52.1 +4.3 −5.1 keV, and the typical νF ν peak energy E peak is also becoming softer, moving from ∼ 1 MeV to 450 keV. The spectral indices of the two power-laws below and above the low-energy break are distributed around the values of -0.82 and -1.72, which are similar to the predictions for synchrotron emission in marginally fast-cooling regime (i.e. −2/3 and −3/2). This is consistent with what has been found in GRB 211211A [20]. We notice, however, that in all spectra the highenergy power-law above E peak is characterised by a much softer index (with a mean value of −4.10 ± 0.24) with respect to the value of ∼ -2.5 typically found in Fermi GRBs. This suggests that the spectral data might require a cut-off at high energy, although further investigations are needed to support this. From 19.5 s until 40.5 s (the last time interval analysed), all the break energies are found to be below 20 keV, close to the GBM low energy threshold. In the same time intervals, the peak energy E peak decreases from 682.4 +3.2 −6.1 to 123.1 +5.4 −4.9 keV, and the index of the power-law below the peak energy is fully consistent (mean value of −1.45 ± 0.06) with the synchrotron predicted value of -1.5.

NTT -Afterglow discovery
Following the refinement of the IPN error box to an area of 30 arcmin 2 [72], we obtained observations of the field of GRB 230307A with the ULTRACAM instrument [76], mounted on the 3.5m New Technology Telescope (NTT) at La Silla, Chile. The instrument obtains images in 3 simultaneous bands, and is optimised for short exposure, low dead-time observations [76]. We obtained 10×20 s exposures in two pointings in each of the Super SDSS u, g and r, bands (where the Super SDSS bands match the wavelength range of the traditional SDSS filters but with a higher throughput; Dhillon et al. 77). The observations began at 01:53:21 UT on 2023-03-09, approximately 34 hr after the GRB. The images were reduced via the HIPERCAM pipeline [77] using bias and flat frames taken on the same night. Visual inspection of the images compared to those obtained with the Legacy Survey [7] revealed a new source coincident with an X-ray source identified via Swift/XRT observations [5], and we identified it as the likely optical afterglow of GRB 230307A [8]. The best available optical position of this source (ultimately measured from our JWST observations, see below) is RA(J2000) = 04:03:26.02, Dec(J2000) = −75:22:42.76, with an uncertainty of 0.05 arcseconds in each axis. The IPN error box and the footprint of the ULTRACAM observations are shown in Figure 5.
This identification was subsequently confirmed via observations from a number of additional observatories, including [9][10][11][78][79][80]. We acquired two further epochs of observations with ULTRACAM on the following nights with 10 × 20s exposures in the Super SDSS u, g and i, bands. Aperture photometry of the source is reported in Table 1, and is reported relative to the Legacy survey for the gri bands, and to SkyMapper for the u-band.

TESS
The prompt and afterglow emission of GRB 230307A was detected by TESS, which observed the field continuously from 3 days before the Fermi trigger to 3 days after at a cadence of 200 s [12]. A reference image was subtracted from the observations to obtain GRB-only flux over this period. The measured flux in the broad TESS filter (600nm -1000nm) is corrected for Galactic extinction and converted to the I c band assuming a power-law spectrum with F ∝ ν −0.8 . We then bin the light curve logarithmically, taking the mean flux of the observations in each bin and converting to AB magnitudes. A systematic error of 0.1 magnitudes was added in quadrature to the measured statistical errors to account for the uncertainties in the data processing. These data are presented in Table 1.

Swift/UVOT
The Swift Ultra-violet/Optical Telescope [UVOT; 81] began observing the field of GRB 230307A ∼ 84.6 ks after the Fermi/GBM trigger [1]. The source counts were extracted using a source region of 5 arcsec radius. Background counts were extracted using a circular region of 20 arcsec radius located in a source-free part of the sky. The count rates were obtained from the image lists using the Swift tool uvotsource. A faint catalogued unrelated source also falls within the 5 arcsec radius, this will affect the photometry, particularly at late times. We, therefore, requested a deep template image in white in order to estimate the level of contamination. We extracted the count rate in the template image using the same 5 arcsec radii aperture. This was subtracted from the source count rates to obtain the afterglow count rates. The afterglow count rates were converted to magnitudes using the UVOT photometric zero points [82,83].

Gemini
We obtained three epochs of K-band observations using the FLAMINGOS-2 instrument on the Gemini-South telescope. These observations were reduced through the DRAGONS pipeline to produce dark and sky-subtracted and flat-fielded images [84]. At the location of the optical counterpart to GRB 230307A, we identify a relatively bright K-band source in the first and second epochs, with only a upper limit in epoch 3. We report our photometry, performed relative to secondary standards in the VISTA hemisphere survey [85], in Table 1.

VLT imaging
We carried out observations of the GRB 230307A field with the 8.2-m VLT telescopes located in Cerro Paranal, Chile. The observations were obtained with the FORS2 camera (mounted on the Unit Telescope 1, UT1, ANTU) in B, R, I and z bands at multiple epochs, and with the HAWK-I instrument (mounted on the Unit Telescope 4, UT4, Yepun) in the K band at one epoch. All images were reduced using the standard ESO (European Southern Observatory) Reflex pipeline [86]. The source was detected in the FORS2 z-band image at ∼6.4 days after the Fermi/GBM detection. A single rband observation of the GRB 230307A was also executed with the 2.6m VLT Survey Telescope (VST) after 2.37 days from the GRB discovery. In later observations the source was not detected (see bottom right panels of Fig. 5) and the upper limit values at 3σ level are reported in Table 1.

VLT spectroscopy
To attempt to measure the redshift of GRB 230307A and the nearby candidate host galaxies, we obtained spectroscopy with the VLT utilising both the X-shooter and MUSE instruments, mounted respectively on the Unit Telescope 2 (UT2, Kueyen) and on UT4 (Yepun).
X-shooter spectroscopy, covering the wavelength range 3000-22000Å was undertaken on 2023-03-15. Observations were taken at a fixed position angle, with the slit centred on a nearby bright star. X-shooter data have been reduced with standard esorex recipes. Given that only two of the four nod exposures were covering the GRB position, resulting in a total exposure time of 2400s on-source, we have reduced each single exposure using the stare mode data reduction. Then, we have stacked the two 2D frames covering the GRB position using dedicated post-processing tools developed in a python framework [87]. We further obtained observations with the MUSE integral field unit on 2023-03-23. The MUSE observations cover multiple galaxies within the field, as well as the GRB position, and cover the wavelength range 4750-9350Å. MUSE data have been reduced using standard esorex recipes embedded within a single python script that performs the entire data reduction procedure. Later, the resulting datacube has been corrected for sky emission residuals using ZAP [88]. The MUSE observations reveal the redshifts for a large number of galaxies in the field, including a prominent spiral G1 at z = 0.0646 [see also 80] and a group of galaxies, G2, G3 and G4 at z = 0.263.

X-ray afterglow
Swift began tiled observations of the IPN localisation region with its X-ray Telescope [XRT; 89] at 12:56:42 on 8 Mar 2023 2 [90]. XRT made the first reported detection of the afterglow (initially identified as 'Source 2') with a count rate of 0.019±0.004 cts −1 [91] and later confirmed it to be fading with a temporal power-law index of 1.1 +0.6 −0.5 [92]. XRT data were downloaded from the UK Swift Science Data Centre [UKSSDC; 93,94].
We further obtained observations with the Chandra X-ray observatory (programme ID 402458: PI Fong/Gompertz). A total of 50.26 ks (49.67 ks of effective exposure) of data were obtained in three visits between 31 March 2023 and 2 April 2023. The source was placed at the default aim point on the S3 chip of the ACIS detector. At the location of the optical and X-ray afterglow of GRB 230307A, we detect a total of 12 counts, with an expected background of ∼ 1, corresponding to a detection of the afterglow at > 5σ based on the photon statistics of [95]. To obtain fluxes, we performed a joint spectral fit of the Chandra and Swift/XRT data. The best fitting spectrum, adopting uniform priors on all parameters, is a power law with a photon index of Γ = 2.50 +0. 30 −0.29 when fitting with a Galactic N H = 1.26 × 10 21 cm −2 [96] and zero intrinsic absorption (neither XRT nor Chandra spectra have sufficient signal to noise to constrain any intrinsic absorption component). The resultant flux in the 0.3 -10 keV band is F X (1.7 d) = 4.91 +0. 89 −0.79 × 10 −13 erg cm −2 s −1 during the XRT observation and F X (24.8 d) = 1.19 +0. 87 −0.62 ×10 −14 erg cm −2 s −1 during the Chandra observation. Due to the low count number, the Chandra flux posterior support extends to considerably below the reported median, with the 5th percentile being as low as F X,5th = 3 × 10 −15 erg cm −2 s −1 . If a uniform-in-the-logarithm prior on the flux were adopted, this would extend to even lower values. Chandra and XRT fluxes are converted to 1 keV flux densities using the best fit spectrum ( Table 2).

ATCA
Following the identification of the optical afterglow [97], we requested target of opportunity observations of GRB 230307A (proposal identification CX529) with the Australia Telescope Compact Array (ATCA) to search for a radio counterpart. These data were processed using Miriad [98], which is the native reduction software package for ATCA data using standard techniques. Flux and bandpass calibration were performed using PKS 1934-638, with phase calibration using interleaved observations of 0454-810.
The first observation took place on 2023-03-12 at 4.46 d post-burst, which was conducted using the 4 cm dual receiver with frequencies centered at 5.5 and 9 GHz, each with a 2 GHz bandwidth. The array was in the 750C configuration 3 with a maximum baseline of 6 km. A radio source was detected at the position of the optical afterglow at 9 GHz with a flux density of 92 ± 22µJy but went undetected at 5.5 GHz (3σ upper limit of 84µJy). A further two follow-up observations were also obtained swapping between the 4 cm and 15 mm dual receivers (the latter with central frequencies of 16.7 and 21.2 GHz, each with a 2 GHz bandwidth). During our second epoch at 10.66 d we detected the radio counterpart again, having become detectable at 5.5 GHz with marginal fading at 9 GHz. By the third epoch, the radio afterglow had faded below detectability. We did not detect the radio transient at 16.7 or 21.2 GHz in either epoch. All ATCA flux densities are listed in Table 2.

MeerKAT
We were awarded time to observe the position of GRB 230307A with the MeerKAT radio telescope via a successful Director's Discretionary Time proposal (PI: Rhodes, DDT-20230313-LR-01). The MeerKAT radio telescope is a 64-dish interferometer based in the Karoo Desert, Northern Cape, South Africa [99]. Each dish is 12 m in diameter and the longest baseline is ∼8 km allowing for an angular resolution of ∼ 7 arccsec and a field of view of 1 deg 2 . The observations we were awarded were made at both L and S-band. GRB 230307A was observed over three separate epochs between seven and 41 days post-burst. The first two observations were made at both L and S4-band (the highest frequency of the five S-band sub-bands), centred at 1.28 and 3.06 GHz with a bandwidth of 0.856 and 0.875 GHz, respectively. Each observation spent two hours at L-band and 20 minutes at S4-band. The final observation was made only at S4-band with one hour on target. Please see the paper by MPIfR for further details on the new MeerKAT S-band receiver.
Each observation was processed using oxkat, a series of semi-automated Python scripts designed specifically to deal with MeerKAT imaging data [100]. The scripts average the data and perform flagging on the calibrators from which delay, bandpass and gain corrections are calculated and then applied to the target. The sources J0408-6545 and J0252-7104 were used at the flux and complex gain calibrator, respectively. Flagging and imaging of the target field are performed. We also perform a single round of phase-only self-calibration. We do not detect a radio counterpart in any epoch in either band. The rms noise in the field was measured using an empty region of the sky and used to calculate 3σ upper limits which are given in Table 2. At the first epoch, we obtained observations in the F070W, F115W, F150W, F277W, F356W and F444W filters of NIRCam [101], as well as a prism spectrum with NIRSpec [102]. In the second epoch we obtained NIRCam observations in F115W, F150W, F277W and F444W and a further NIRSpec prism observation. However, in the second epoch the prism observation is contaminated by light from the diffraction spike of a nearby star and is of limited use, in particular at the blue end of the spectrum. We therefore use only light redwad of 1.8 microns. However, even here we should be cautious in interpreting the overall spectral shape. However, the feature at 2.15 microns is visible in both the 29 and 61 day spectra.

JWST observations
We reprocessed and re-drizzled the NIRCam data products to remove 1/f striping and aid point spread function recovery, with the final images having plate scales of 0.02 arcsec/pixel (blue channel) and 0.04 arcsec/pixel (red channel).
In the NIRCam imaging we detect a source at the location of the optical counterpart of GRB 230307A. This source is weakly detected in all three bluer filters (F070W, F115W and F150W), but is at high signal-to-noise ratio in the redder channels (see Figure 2). The source is compact. We also identify a second source offset (H1) approximately 0.3 arcseconds from the burst location. This source is also weakly, or non-detected in the bluer bands, and is brightest in the F277W filter.
Because of the proximity of the nearby star and a contribution from diffraction spikes close to the afterglow position we model point spread functions for the appropriate bands using WebbPSF [103], and then scale and subtract these from star position. Photometry is measured in small (0.05 arcsec (blue) and 0.1 arcsec (red)) apertures and then corrected using tabulated encircled energy corrections. In addition to the direct photometry of the NIRCam images we also report a a K-band point based on folding the NIRSpec spectrum (see below), through a 2MASS, Ks filter. This both provides a better broadband SED and a direct comparison with ground based K-band observations. Details of photometric measurements are shown in Table 1.
For NIRSpec, we utilise the available archive-processed level 3 two-dimensional spectrum ( Figure 3). In this spectrum we clearly identify the trace of the optical counterpart, which appears effectively undetected until 2 microns and then rises rapidly. We also identify two likely emission lines which are offset from the burst position. These are consistent with the identification with Hα and [O iii] (4959/5007) at a redshift of z = 3.87. Both of these lines lie within the F277W filter in NIRCam and support the identification of the nearby source as the origin of these lines.
We extract the spectrum in two small (2 pixel) apertures. One of these is centred on the transient position, while the second is centred on the location of the emission lines. Since the offset between these two locations is only ∼ 0.2 arcseconds there is naturally some contamination of each spectrum with light from both sources, but this is minimised by the use of small extraction apertures. The resulting 1D spectra are shown in Figure 3. The counterpart is very red, with a sharp break at 2 microns and an apparent emission feature at 2.15 microns. The spectrum then continues to rise to a possible second feature (or a change in the associated spectral slope) at around 4.5 microns. 20 3 Supplementary Information 4 GRB 230307A in context 4.1 Prompt emission GRB 230307A is an exceptionally bright GRB. It has the second highest fluence of any GRB observed in more than 50 years of GRB observations [52]. While it remains a factor of 50 less fluent than GRB 221009A, it is still a factor ∼ 2 brighter than GRB 130427A, the third brightest burst. Bursts with these extreme fluences are rare. In Figure 6, we plot the distribution of observed fluence for Fermi/GBM detected bursts. At the brighter end, the slope of the distribution is consistent with the expected −3/2 slope for a uniform distribution of sources. The extrapolation of this relation suggests that bursts like GRB 230307A should occur once every several decades. Notably, three bursts well above the extrapolation (GRB 130427A, GRB 230307A, GRB 221009A) may indicate that bright bursts arise more frequently than expected. However, observationally it is clear that GRB 230307A is, at least, a once-per-decade event.
The prompt light curve of GRB 230307A ( Figure 1) shows two distinct emission features: an initial episode of hard emission from the trigger until ≈ 18 s, then a softer episode from ≈ 19 s onwards. These distinct episodes of hard and soft emission are strongly reminiscent of the long-duration merger GRB 211211A, but the initial pulse complex is ∼ 50 per cent longer in GRB 230307A when compared to the ∼ 12 s duration seen in GRB 211211A [17]. The relative durations of the initial pulse complex in the two GRBs bear a striking resemblance to their relative time-averaged peak energies [936 ± 3 keV vs 647 ± 8 keV; 2, 104]. In GRB 211211A, substantial spectral evolution was seen to drive the light curve, and the underlying radiation mechanism was identified as fast-cooling synchrotron emission [20]. The coherent development of the hardness ratio ( Figure 1, lower) indicated similar spectral evolution in GRB 230307A, which the spectral analysis confirmed. Indeed, as described in Section 2.1.1, the time-resolved spectral analysis of the prompt emission revealed the presence of two spectral breaks in the GBM band, E break and E peak , coherently becoming softer from 7.5 s up to 19.5 s. Also, in this case, the spectral indices indicate synchrotron emission in the marginally fast-cooling regime. From 19.5 s onwards (approximately when the softer and dimmer emission episode starts), the low-energy break E break is continuously approaching the lower limit of the GBM band (8 keV), presumably crossing it to enter the X-ray regime. Unfortunately, the lack of simultaneous observations in X-rays with another telescope, e.g. Swift/XRT, prevents us from fully tracing the evolution of the spectral break down to X-rays at later times, as was done for GRB 211211A.
The time-averaged Fermi/GBM spectrum of GRB 230307A across the T 90 interval is best fit with a cutoff power-law with α = 1.07 ± 0.01 and cutoff energy 936 ± 3 keV [2]. From this, we calculate a hardness ratio (the ratio of the 50 -300 keV photon flux to the 10 -50 keV photon flux) of 0.88 +0.01 −0.02 . This is higher than the value for 211211A (0.57) but comfortably within the 1σ distribution of hardness ratios for canonical long GRBs (i.e. with T 90 > 2 s) in the Fermi catalogue, which we calculate to be 0.66 +0. 51 −0.29 from the data in von Kienlin et al. [105]. Like GRB 211211A before it, GRB 230307A appears to have 'typical' long GRB properties in terms of its time-averaged hardness ratio and its T 90 . This strengthens the case for a significant number [17,18] of long-duration GRBs having been mistakenly identified as stellar collapse events. However, in some ways, GRB 230307A differs significantly from several of the other brightest GRBs. For example, the afterglow was relatively faint, while the burst was very bright. In Figure 7, we plot the prompt fluence in the 15-150 keV band against the X-ray afterglow brightness at 11 hours [updated from 106,107]. The general trend between the afterglow brightness and fluence is seen; the best-fit slope to this relation is approximately one. So, while there is substantial scatter, there is a direct proportionality between the fluence and the afterglow brightness. Notably, while the afterglow and prompt emission of GRB 221009A were exceptionally bright (after correcting for the heavy foreground extinction), they were in keeping with this relatively broad relationship. GRB 230307A is different. Here we extrapolate the X-ray flux to 11 hours based on the measured X-ray flux at ∼ 1 day and the decay slope. We also re-calculate the GRB 230307A fluence in the relevant 15-150 keV energy band for comparison to Swift/BAT. This burst is a notable outlier in the relation, with a faint X-ray flux for its extraordinary prompt brightness. The afterglow brightness depends both on the energy of the burst and the density of the interstellar medium; it is, therefore, possible that the location in this fluence -afterglow brightness plane is indicative of a lowdensity medium, which would be consistent with expectations for such a large GRBhost offset.
It is also of interest that another burst in a similar location is GRB 211211A. This long burst has a clear signature of kilonova emission within its light curve. If GRB 230307A is a similar event, faint afterglows (relative to the prompt emission) may be an effective route for disentangling mergers from collapsars.
To further compare the ratio between the X-ray brightness and the γ-ray fluence, we retrieve the X-ray light curve of all Swift-detected GRBs from the Swift Burst Analyser [108] and limit the sample to 985 long GRBs and 55 short GRBs with at least two XRT detections and measured BAT fluence. The fluences are taken from the Swift/BAT Gamma-Ray Burst Catalog 4 [109] and represent the measurements from 15 to 150 keV integrated over the total burst duration. We add to this sample the GRBs 170817A [off-axis GRB; 110] and 221009A [brightest GRB detected to date; 52]. For the former, we retrieve the X-ray light curve from Hajela et al. [111] and use a γray fluence of 2.4 × 10 −7 erg cm −2 [110]. For the latter, we take the X-ray light curve from the Swift Burst Analyser and assume a fluence of 0.007 erg cm −2 (corrected from the 1-10000 keV fluence in [52] to the 15-150 keV band). Following [112], we resample the X-ray light curves and normalise them by the γ-ray fluence on a grid defined by the observed F X /Fluence ratios and the time-span probed by the data. If no data are available at a specific time of the grid, we linearly interpolate between adjacent observations but do not extrapolate any data. Hence the paucity of observations at later times reflects the last time at which sources were detected by the Swift/XRT.
Short and long GRBs occupy the same part of the F X /Fluence vs time parameter space ( Figure 1). In contrast, GRB 230307A has an unprecedentedly low F X /Fluence ratio that is almost 10-fold lower than the faintest GRBs at the same time. To emphasise the uniqueness of GRB 230307A, we also show in the same figure the Swift/BAT-detected GRBs 050925, 051105A, 070209, 070810B, 100628A, 130313A, 170112A that evaded detection with Swift/XRT. The limits on their F X /Fluence ratio (shown by downward pointing triangles in that figure) are consistent with the observed range of F X /Fluence ratios, ruling out a selection bias against GRBs with lower than usual F X /Fluence. Intriguingly, GRBs 080503, 191019A and 211211A had markedly low F X /Fluence ratios during the shallow decline phase of their X-ray light curves. Furthermore, GRB 211211A reached a value of 1.2 × 10 −9 s −1 at 120 ks, comparable to GRB 230307A.

Counterpart Evolution
Although the afterglow of GRB 230307A was promptly detected thanks to TESS, this data was not available to the community for several days. Further follow-up was, therefore, much slower, and the counterpart was not discovered until the localisation was narrowed down to several sq. arcminutes, approximately 24 hours after the burst. The result is that the counterpart is poorly sampled (particularly in colour) during the early phases, while later observations suffer from typically modest signal-to-noise.
The TESS observations detected a relatively bright (though not exceptional given the fluence of the burst) outburst, coincident with the prompt emission, likely peaking at I < 15 [12]. The afterglow was much fainter, apparently no brighter than I = 18 in the minutes to hours after the burst was detected. It was relatively flat during this period, with a power-law through the first to last TESS observations decaying as F (t) ∝ t −0.2 . The TESS and ground-based observations can be consistently modelled with a forward shock afterglow + kilonova (see Section 6.1).
There are no simultaneous colours at the time of the first ground-based afterglow detections (1.4 days), although extrapolation of the r-band detection with ULTRA-CAM to the WHITE detection with the Swift/UVOT suggests a relatively red colour (WHITE-r = 1.6 ± 0.4) . However, such an interpretation is difficult due both to the large photometric errors and the width of the WHITE filter on the Swift/UVOT.
Optical observations obtained multiple colours at an epoch ∼ 2.4 days post-burst. These show the afterglow to have a blue colour with g = 22.35 ± 0.26, i = 21.68 ± 0.09 and z = 21.8 [80]. This is consistent with GRB afterglows in general(i.e F ν ∝ ν −β gives β ≈ 1). Observations in the near infrared (NIR) were not undertaken until ∼ 10.4 days post-burst. However, these reveal a relatively bright K-band source. The inferred i-K(AB) > 2.9 at this epoch is very red. Interpreted as a change in the spectral slope, it is β ≈ 2.5. The K-band light hence appears to be in significant excess with respect to the afterglow expectations based both on optical data and on the X-ray light curve.
It is relevant to consider if such an excess could arise via extinction. However, this is not straightforward to explain. For a generic β = 1 slope we expect i−K(AB) ≈ 1.1. At z = 0, to obtain i − K(AB) = 2.9 would correspond to a foreground extinction of A V ≈ 4. However, this would also predict g − i ≈ 3, which is entirely inconsistent with the earlier observations. This problem becomes more acute for higher redshifts, where the bluer bands probe increasingly into the UV.
The IR excess becomes extremely prominent by the time of the JWST observations. At 28.5 days, the source is detected in all bands but is very faint in the NIRCam blue channel (F070W, F115W, F150W) and rises rapidly (in F ν ) through the redder bands (F277W, F356W, F444W). Expressed as a power-law, this is β ≈ 3.1 in the 2-5 micron region, and β ≈ 1 between 0.7-1.5 microns. This does not match the expectations for any plausible spectral break in a GRB afterglow or any plausible extinction (where one would expect the slope to steepen towards the blue). This strongly implies that the red excess seen in the K-band at ten days and with JWST at 28.5 days is some additional component. Indeed, in the JWST observations, the other component, beginning at around 2 microns, is very clearly visible in both photometry and spectroscopy.
This component evolves exceptionally rapidly. In the K-band, the inferred decay rate from 11.5 to 28.5 days is ∼ t −3.5 expressed as a power-law or ∼ 0.25 mag per day, if exponential. This is much faster than observed in GRB afterglows or supernovae. It is, however, consistent with the expectations for kilonovae. As shown in Figure 4, the overall evolution shows substantial similarity with AT2017gfo. To constrain the temporal and spectral evolution within a plausible physical model more accurately, we fit the multi-band photometry with afterglow and kilonova models. The outputs of these models are described in detail in section 6.

Identification of the host galaxy
Deep optical imaging of the field identifies several relatively bright galaxies in the vicinity of the sky position of GRB 230307A. Our preferred host galaxy is the brightest of these, which we denote as G1. It lies at z = 0.065 and is offset 30 arcseconds (40 kpc in projection) from the location of the afterglow. Following the method of [113] this galaxy has a probability of chance alignment of P chance ∼ 0.09 (see also [80]). Although this is not extremely low, and so is only suggestive of a connection to the transient, we note that i) the luminosity of the late time counterpart at this redshift is very similar to AT2017gfo and ii) the spectral feature seen at 2.1 microns in AT2017gfo matches with the emission feature seen in the JWST spectroscopy of GRB 230307A. This is a broad line, but assuming they have the same physical origin, they fix the redshift to the range 0.04 < z < 0.08. G1 is the only galaxy within this range in the field. The physical properties of this galaxy are outlined in section 4.4.
Our MUSE observations provide redshifts for this galaxy and several others, also identifying a small group of galaxies (G2, G3, G4) at a common redshift of z = 0.263. All of these galaxies have P chance values substantially greater than our preferred host. Furthermore, because of the larger redshifts, the implied offsets from GRB 230307A are 100kpc. This is larger than seen for any short GRB with a firmly identified host. We, therefore, disfavour these as plausible host galaxies for GRB 230307A.
Deep JWST observations reveal no evidence of a directly underlying host galaxy for GRB 230307A, as would be expected if it had a collapsar origin. In particular, at late times, the faint source at the counterpart's location is consistent with a pointsource (i.e. a subtraction of the PSF constructed by WebbPSF yields no significant residuals). However, we identify a faint galaxy, undetected in the blue and with F277W =27.9 ± 0.1, offset only 0.3 arcseconds from the burst position. We designate this galaxy H1.
Our NIRSpec observations provide a redshift of z = 3.87 for H1 based on the detection of [O III] (5007) and Hα. At this redshift, the offset is only ∼ 1.3 kpc.
Although many z ∼ 4 galaxies are extremely compact [114], it seems likely that some stellar population from this galaxy does extend under the burst position, and there may be marginal evidence for extension in this direction in the F444W image. However, this region is neither UV-bright nor an emission line region where one may expect to observe massive stars.
The galaxy photometry, performed in 0.1-arcsecond apertures and subsequently corrected for encircled energy assuming point-source curves is F070W>29.0, F115W=28.4 ± 0.3, F150W=28.6 ± 0.4, F277W=27.9 ± 0.1, F444W=28.3 ± 0.1, and the galaxy is only robustly detected in the redder bands (see Figure 2). We note that because of the proximity of the afterglow, we use a smaller aperture than may be optimal, although the galaxy is also compact.
We can estimate the probability of chance alignment of this source with the GRB position via various routes. In principle, one can use number counts of galaxies on the sky in the multiple bands. These have recently been updated based on the first observations with JWST to provide number counts in appropriate bands [115]. We find that P chance , following the approach of [113] to be in the range ∼ 3-6 % for F277W and F444W (with no bound in the filters where the galaxy is undetected). Alternatively, we also estimate the probability directly from the data. We extract sources within the field via Source Extractor to create a mask of objects within the field. In the brightest detection (F277W) approximately 5% of the image is covered with objects of equal or brighter magnitude to H1, and we note that the burst position is not contained within this mask. This suggests that in this particular field, P chance > 5%.
The absolute magnitude of H1 is M i ∼ −17.7, and the Hα star formation rate is approximately 1 M yr −1 . The half-light radius of the galaxy is approximately 0.1 arcseconds (700 pc). Although limited information is available, these values are generally consistent with those of the long GRB population. The burst offset from its host galaxy is ∼ 2.5 half-light radii. This is large but within the range seen for long-GRBs [116].
In our X-shooter and MUSE observations there is no trace visible in 1D or 2D extractions at the source position, although a weak continuum is seen in the X-shooter spectrum when heavily binned. This is consistent with its faint magnitude at the time of the observations. At the location of Lyα at z = 3.87 we place limits of F < 2.5 × 10 −17 erg s −1 cm −2 assuming an unresolved line.
We also examined both spectra for any emission lines at other redshifts. This is worthwhile given the strong emission lines often seen in long GRB hosts [117], which may make emission line redshifts possible, even if the host itself is undetected. However, despite deep observations, there are no visible emission lines consistent with no directly underlying host galaxy, consistent with a compact object merger, but not a collapsar.
Unsurprisingly, there are also numerous faint galaxies in the JWST images. However, all of these have large P chance values, and we do not consider them plausible host galaxies.
Taken a face value, the probability of chance alignment for G1 (our preferred host) and H1 (z = 3.87) is similar. However, the luminosity, lightcurve evolution and spectroscopic feature at the redshift of G1 offer strong support for it as the host galaxy of GRB 230307A. Furthermore, there is no straightforward, reasonably viable physical model that could explain the burst's extreme properties at z = 3.87. This scenario would require extreme energetics, exceptionally rapid evolution and yields unphysical outcomes in standard GRB or supernovae scenarios. We outline this in detail in section 7.1.

Host galaxy properties
To better understand the properties of G1, the likely GRB host galaxy, we performed a fit to both the MUSE spectrum and photometric measurements from the far-UV to the mid-IR. For the photometric measurements, we retrieved science-ready coadded images from the Galaxy Evolution Explorer (GALEX) general release 6/7 [118], DESI Legacy Imaging Surveys (LS; [119]) data release 9, and re-processed WISE images [120] from the unWISE archive [121] 5 . The unWISE images are based on the public WISE data and include images from the ongoing NEOWISE-Reactivation mission R7 [122,123]. We measured the brightness of the galaxy G1 using the Lambda Adaptive Multi-Band Deblending Algorithm in R (LAMBDAR [124]) and the methods described in [125]. We augment the SED with Swift/UVOT photometry in the u band and our 6-band JWST/NIRCAM photometry. The photometry on the UVOT images was done with uvotsource in HEASoft and an aperture encircling the entire galaxy. For JWST photometry, we used a 6-arcsec circular aperture, which allows us to gather all the observed light observed in JWST filters from the host galaxy. All measurements are summarised in Table 7.4.
To derive the main physical properties of the host galaxy, such as its stellar mass, we employ two separate methodologies based on the photometric and spectroscopic data available for the host, and finally compare the results to assess the robustness of our conclusions. We first fit the multi-wavelength (0.1-4.4 µm) dataset using the prospector python package [126], which allows us to model the host galaxy spectrum starting from its main constituents, namely a set of stellar population base spectra, built from the Flexible Stellar Population Synthesis (FSPS) package [127], and combined with a specific star-formation history (SFH) model. Moreover, we have also considered a fixed attenuation model based on the Calzetti [128] attenuation curve, and an additional nebular model originating from the gas component, which is built using the Cloudy photo-ionization code [129], and considering the FSPS stellar population as ionising sources. We have adopted a parametric SFH model, which is described by a delayed-exponential model where the star-formation rate varies as a function of time t = t age − t lt , with t lt being the lookback time [126], as SFR ∝ (t/τ ) exp(−t/τ ), with τ being the e-folding time. We finally have used the dynesty Speagle [130] ensemble sampler to reconstruct the posterior distribution.
The results of the prospector analysis are shown in Fig. 8. We obtain a mass value of the living stars of M * = 2.37(+0.24, −0.35) × 10 9 M yr −1 . The mass of all stars ever formed is 0.20(+0.02, −0.04) dex larger. The light-weighted stellar age resulting from the fit is 1.13(+1.49, −0.36) Gyr.
An alternative to parametric SED fits is to use synthetic stellar population SEDs as templates and combine them to fit the galactic spectra (the underlying assumption being instantaneous star formations rather than continuous functions of time). We can use the spectral synthesis from the BPASSv2.2.2 [131,132] binary populations and create templates with hoki [133] that are compatible with the ppxf fitting package [134], as described in [135]. Because SED fitting has a high level of degeneracies (see [134]), at first we do not fit all 13 BPASS metallicities at once with ppxf, as this can result in unphysical results (see discussion in Stevance et al. 135); instead we fit the metallicities individually to find which ones result in the best fits on their own. We find that a low Z (0.001) population and solar metallicity population (Z=0.014) result in decent fits, but the low metallicity population fails to predict a young stellar component that is seen in the images, whilst the solar metallicity fit fails to accurately match the Hβ and neighbouring absorption features in the blue part of the spectrum. So we then fit the galaxy simultaneously with Z=0.001 and Z=0.014 templates, and retrieve a good fit shown in Figure 9 alongside the recovered SFH.
We find evidence of three main stellar populations: >95% of the mass is found in lower metallicity (Z=0.001) stars with ages ranging from a few Gyr to 10 Gyr, with a peak of star formation around 5 Gyr; >4.7% of the mass originates from a solar metallicity population (Z=0.014) that formed around 400 Myr ago; finally a small fraction (<0.05%) of the stellar mass in the host originates from the star-forming regions with ages a few Myrs.
The details of the age distributions and exact metallicity values can be model dependent so we also fit the integrated galaxy spectrum with the single stellar population synthesis code STARLIGHT [136], which uses stellar populations based on 25 different ages and six metallicity values [137], and a Chabrier IMF [138]. The SFH retrieved by this method is more complex and would require odd configurations (including some high metallicties at old ages and low metallicities around 100 Myr, which is counter-intuitive, unless inflow from pristine gas will trigger a burst of SFR), but it also finds that overall the galaxy is dominated by an old population with lower metallicity and has a younger component at higher metallicity. In Figure 9 we show a comparison of the STARLIGHT and BPASS fits in the bottom left panel and see that they are very similar, despite STARLIGHT containing 6 different metallicities and assuming solely single star populations. This highlights the level of degeneracy we face when performing galaxy SED fits. We leave further comparisons to a follow-up study dedicated to the host and the progenitor populations of GRB 230307A, where we will also present detailed, specially resolved, fits to the datacube including its kinematics.
For now we use the BPASS integrated fits to infer the stellar mass and the star formation rate of the host of GRB 230307A, as the fit and SFH is more convincing that the one obtained with STARLIGHT. We find that there are currently M * = 1.65 × 10 9 M of living stars (corresponding to 3.1 ×10 9 M at ZAMS) in G1. Using the nebular component retrieved from subtracting the fit of the stellar component to the observed data, we can also estimate the star formation rate and metallicity. From the Hα feature we estimate that the SFR is 5.47±0.30 × 10 −1 M yr −1 using the Kennicutt formulation [139], and using the N2 index, in the CALIFA formulation [140], we infer an oxygen abundance of 8.20 ± 0.16 (12 + log(O/H)).
There are qualitative similarities between the host of GRB 230307A and NGC 4993 [135], the host of the first confirmed kilonova (they are both dominated by an older stellar populations and include a younger more metal rich component), but there 27 are some key differences: NGC 4993 was a lenticular galaxy without a clear young component, whereas the host of GRB 230307A shows clear spiral arms and star forming regions. Another major difference is that the metallicity of the old population in this galaxy is 10 times lower than that of NGC 4993 (Z=0.001 compared to Z=0.010), which will influence the stellar evolution of potential progenitors. Finally, NGC 4993 had a large stellar (and presumably dark halo) mass M * ≈ 10 11 M [135], a factor of > 50 larger than the host of GRB 230307A.
The location of GRB 230307A relative to its host galaxy is consistent with these properties. In particular, the low mass of the galaxy suggests a modest gravitational potential such that binaries with velocities of a few hundred km s −1 can readily escape. The large offset also suggests that the binary is formed from the older stellar population.

Properties of the brightest GRBs
GRB 230307A is the second brightest 6 burst observed in over 50 years of observations [52]. If it arises from a compact object merger, this implies that such bright bursts can be created in mergers. Indeed, such a picture appears likely based on GRB 211211A [17][18][19], the sixth brightest burst. Of the ten brightest bursts observed by the Fermi/GBM, and subsequently localised at the arcsecond level, three have apparently secure associations with supernovae (GRB 130427A, GRB 171010A, GRB 190114C), and two (GRB 211211A, GRB 230307A) are associated with kilonovae, and hence mergers. Of the remaining five, one lies at z = 1.4 and has energetics which suggest a collapsar; three have no redshift information, although one of these (GRB 160821A) lies in proximity to several galaxies at z = 0.19; and one is GRB 221009A whose associated with a supernova remains unclear [97,141,142], although recent observations suggest a collapsar with an associated supernova is most likely [143]. Within this very bright population, collapsars are likely as common as mergers.

Event rates
One key question of interest is the likely event rate for such merger GRBs. A simple estimate of the event rate associated to a single event is given by Here Ω reflects the fraction of the sky covered by the detection mission, t the effective mission duration (accounting for the duty cycle) and V max the maximum co-moving volume within which a burst with the same properties could be identified. For GRB 230307A, Ω = 0.65 (average for the Fermi/GBM) and t ≈ 15 years. V max is more complicated: as shown in Figure 6, the fluence distribution for GBM bursts extends to ∼ 10 −8 erg cm −2 and is likely complete to around 10 −6 erg cm −2 . Given the extreme brightness of GRB 230307A, it would likely have been recovered to a distance ∼ 50 times greater than its observed distance. If at z = 0.065 the inferred z max = 2.03 6 We use brightest here as an indicator of the total fluence in the prompt emission or V max = 630 Gpc −3 . In this case, the inferred rate of such bursts becomes extremely small, R ≈ 1.6×10 −4 Gpc −3 yr −1 . However, in practice, such bursts would not readily be identified at such redshifts since neither supernova nor kilonova signatures could be observed. A more realistic estimate would correspond to the distance at which associated supernovae can be either identified or ruled out with moderate confidence. In this case z max = 0.5 (also adopted by [144]), V max = 29 Gpc 3 , and R ≈ 3.5 × 10 −3 Gpc −3 yr −1 .
These rate estimates also assume that GRB 230307A is the only merger-GRB to have occurred within the 15-year lifetime of the Fermi/GBM. This is almost certainly not the case. Indeed, GRB 211211A was also identified by Fermi/GBM and has rather similar estimates of the intrinisc rate [5.7 × 10 −3 Gpc −3 yr −1 , 144].
However, even the interpretation of ∼ 2 events is problematic. In particular, the V /V max for GRB 230307A is 0.004, and for GRB 211211A = 0.005 (again assuming z max = 0.5). For a sample average of uniformly distributed sources of comparable energy or luminosity, we expect V /V max ∼ 0.5). That the initial identification of such a population should arise from bursts with such extreme V /V max values is surprising, but may reflect that these bursts are the brightest, which likely encouraged a detailed follow-up. However, it is improbable that they represent the only such bursts observed, and we should expect a much larger population.
To better quantify this, we extend our analysis to the Swift bursts and utilise the fluence of GRB 230307A converted to a 15-150 keV equivalent fluence using the observed spectral parameters. At z = 0.065, E iso (15-150) keV ∼ 7 × 10 51 erg, and for GRB 211211A E iso (15-150) keV = 2 × 10 51 erg. As expected, low energy events dominate the low redshift GRB population. However, at z < 0.5, there are 12 (out of 42) bursts with E iso 10 51 erg. This includes some further supernova-less GRBs, in particular GRB 060614 (E iso = 9 × 10 50 erg), GRB 191019A (E iso = 2.0 × 10 51 erg), and some bursts for which supernova searches have not been reported (e.g. GRB 150727A, GRB 061021, and the 'ultra-long' GRB 130925A). This sets an upper limit on the number of bursts at low redshift, which may be associated with mergers. In practice, selection effects would support a scenario where mergers generate a larger fraction of these bursts. In particular, the afterglows of GRB 230307A and GRB 211211A appear to be faint, despite the bright prompt emission. Such afterglows are difficult to find and may evade detection. In these cases redshifts may only be obtained from host galaxies. The associations may not be obvious if the bursts are offset from host galaxies at moderate redshifts. Such follow-up may occur late after the burst, or optical afterglow non-detections may lead to a lack of optical/IR follow-up because of uncertainty regarding the optical brightness of the event or suggestions it may be optically dark because of host galaxy extinction. Finally, given the afterglow brightness issues, it is possible that the small fraction of bright GRBs without redshift measurements may arise from a similar channel. These observations would imply that between 30-70% bursts at z < 0.5 and E iso 10 51 erg could arise from mergers, although it is likely less. A modest number of events at higher redshift is consistent with the observations, and would alleviate concerns regarding V /V max for GRB 211211A and GRB 230307A.
This fraction is surprisingly high given the strong evidence that long GRBs arise from broad-lined type Ic supernovae and short GRBs from compact object mergers.
However, the dominant contributors to the long-GRB supernova connection occur at low energy, and belong to a population of low luminosity GRBs (LLGRBs) [145]. In a significant number of these, we may observe a energy source in the prompt emission separate from the highly relativistic jet seen in on-axis, energetic bursts. For example, the long-lived, soft nature of some bursts suggests a contribution from shock breakout or cocoon emission. If, for this reason, the luminosity function of collapsar GRBs is steeper at low luminosity than that of merger-GRBs, it is possible that at low luminosity the long GRB population is dominated by collapsars, while at high luminosity the contribution of mergers is significant. Such an interpretation is not without problems, given the star-forming nature of long-GRB hosts and their typically small offsets from their host galaxies. However, it is a logical investigation for future work.

Light curve modelling
In order to shed light onto the properties of the jet and, even more importantly, to separate the contribution of the kilonova from that of the jet afterglow in the UVOIR bands, we modelled the multi-wavelength light curves from radio to X-rays as a superposition of synchrotron emission from the forward shock driven by the jet into the interstellar medium (ISM), following [32,146], and blackbody emission from the photophere of a kilonova, using the simple single-component model of [147].
The forward shock synchrotron emission model has eight parameters, namely the isotropic-equivalent kinetic energy in the jet E K , its initial bulk Lorentz factor Γ 0 , its half-opening angle θ j , the ISM number density n, the fraction ξ N of ISM electrons that undergo diffusive shock acceleration in the forward shock, the fraction e of the shock downstream internal energy that is shared by such electrons, the slope p of the power law dN e /dγ ∝ γ −p that describes the Lorentz factor (as measured in the shock downstream comoving frame) distribution of the accelerated electrons as they leave the acceleration region, and the fraction B of the shock downstream internal energy that is shared by a small-scale, turbulence-driven, random magnetic field. The shock hydrodynamics is computed from energy conservation and accounts for the lateral expansion of the shock [146]. The effective electron energy distribution is computed accounting for the cooling induced by synchrotron and synchrotron-self-Compton emission, including an approximate treatment of the Klein-Nishina suppression of the Thomson cross section [146]. In computing flux densities, the synchrotron surface brightness of the shock is integrated over equal-arrival-time surfaces to account for the effects of relativistic aberration and latitude-dependent retarded times on the spectral shape [148].
The kilonova model [147] assumes spherical ejecta expanding homologously, v = r/t, and featuring a power law density profile ρ(r, t) ∝ t −3 v −δ between a minimum and a maximum velocity, v ej ≤ v ≤ v ej,max . The density normalization is set by the total ejecta mass M ej . In general, the model allows for the ejecta opacity (assumed grey) κ to be piecewise-constant within the profile, but here we assume a uniform opacity across the ejecta for simplicity. The model divides the ejecta into 100 small shells and computes the heating rate and thermalization efficiency within each. This allows for the derivation of the internal energy evolution in each shell and eventually the computation of the photospheric luminosity L KN in the diffusion approximation. The fixed ejecta opacity also allows for the computation of the optical depth and hence for the identification of a photospheric radius, which then sets the effective temperature T KN by the Stefan-Boltzmann law. In our modelling of GRB 230307A, we computed the flux density by simply assuming pure blackbody emission with the given luminosity and effective temperature at each given time. We fixed v max = 0.6c and left M ej , v ej , κ and δ as free parameters.
To carry out the model fitting, we defined an asymmetric Gaussian log-likelihood term for the i-th datapoint, which corresponds to an observation at time t i and in a band whose central frequency is ν i , as where F ν,m (ν, t) is the flux density predicted by the model, F ν,obs,i is the measured flux density, the one-sigma error reflects the potentially asymmetric error bars and we introduced a fractional systematic error contribution f sys , which we take as an additional nuisance parameter, to account for potential inter-calibration uncertainties between different instruments and for the fact that error bars typically only account for statistical uncertainties. For X-ray detections, we fit the integrated flux and the spectral index independently, with an analogous term for each (but with no systematic error contribution for the spectral index). Upper limits were treated simply by setting F ν,obs,i equal to the reported upper limit, σ h,i = F ν,obs,i /10 and σ l,i = 10F ν,obs,i . The final log-likelihood was taken as the sum of these terms. In order to derive a posterior probability density on our 13-dimensional parameter space, we assumed the priors reported in Table 6 and we sampled the posterior with a Markov Chain Monte Carlo approach using the emcee python module [149], which implements the Goodman and Weare [150] affine-invariant ensemple sampler. The medians and 90% credible intervals of the marginalised posteriors on each parameter obtained in this way are reported in Table 6. The posterior is visualised by means of corner plots in Figures 10 (jet afterglow parameters), 11 (kilonova parameters) and 12 (all parameters).
The left-hand panel in Figure 13 shows the observed light curve data (markers) along with the best-fitting model (solid lines). Dashed lines single out the contribution of the kilonova. The right-hand panel in the same figure shows some selected spectra, showing in particular the good agreement of the first JWST epoch with the blackbody plus power law spectrum implied by our model at those times.
While the best-fit model demonstrates a relatively good agreement with most of the measurements, some discrepancies stand out, most prominently with the 61.5 d JWST data and with the 28.5 d Chandra detection. The former is not too surprising, as the assumptions in the kilonova model (in particular that of blackbody photopsheric emission, which is particularly rough in such a nebular phase, and that of constant and uniform grey opacity, due to recombination of at least some species) are expected to break down at such late epochs. The latter is linked to the steepening ('jet break') apparent at around 2 days in the model X-ray light curve, which in turn is mainly driven by the need to not exceed the optical and near-infrared fluxes implied by observations at around one week and beyond. In absence of these constraints, the fit would have accommodated a larger jet half-opening angle, postponing the jet break and hence allowing for a better match with the best-fit Chandra flux. On the other hand, as noted in Methods, this flux is rather uncertain, with the low-end uncertainty possibly extending to fluxes lower by one order of magnitude or more, depending on the adopted prior in the spectral analysis (see Methods). Still, such a discrepancy might indicate the presence of additional X-ray emission that is not accounted for by the model, as has been seen previously in e.g. [151,152].

Spectral analysis modeling
The JWST/NIRSpec spectrum taken on 5 April 2023 exhibits a red continuum component with emission line features. The most distinctive feature is a broad emission line at 2.15 microns (in the rest frame, assuming z = 0.065). This may be a blend (visibly split in Figure 3) and a simultaneous fit of two Gaussians provides measured centroids of 20285±10Å and 22062±10Å. The line widths are both consistent at v FWHM = 19100 kms −1 (0.064c). This 2.1 micron feature is quite similar in strength and width to the 2.07 micron feature in AT2017gfo at 10.5 days after merger [discussed in 37]. The AT2017gfo line also appears to be better fit as a blend of two features rather than a single transition, with line velocities of v FWHM = 38900 kms −1 . While the average line centre is reasonably consistent between the two, the components inferred for AT2017gfo and the kilonova of GRB 230307A are each quite different. Reference [37] finds them at 20590Å and 21350Å and there is no consistent velocity shift that could be applied to match AT2017gfo with our JWST spectrum. Nevertheless, the similarity in their average line centroids, velocities and equivalent widths is striking, as demonstrated in Figure 3.
With a Doppler broadening parameter of 0.1c, it is unlikely that the continuum component is formed as a result of the superposition of emission lines. Because kilonova radiation transfer at such late times is not yet fully understood, here we attempt to model the spectrum with the assumption that the emission consists of blackbody radiation from the photosphere and forbidden emission lines of heavy elements formed outside the photosphere.
If the continuum is described with blackbody radiation, the temperature and photospheric velocity are ≈ 670 K and ≈ 0.08c, respectively. The continuum luminosity is estimated as ∼ 2 × 10 39 erg/s in the NIRSpec band and ∼ 5 × 10 39 erg/s if the blackbody emission extends to much longer wavelengths. Assuming this emission is entirely powered by radioactivity of r-process nuclei, these correspond to an ejecta mass of ∼ 0.03-0.07M [147]. With the ejecta mass and velocity, the opacity is required to be 5 cm 2 /g in order to keep the ejecta optically thick at 30 day. It is worth noting that such a high opacity in the mid-IR indicates that the inner part of the ejecta is lanthanide rich [153][154][155].
Forbidden emission lines in the infrared are expected to arise from fine structure transitions of low-lying energy levels of heavy elements. Most abundant ions are expected to produce the strongest lines. We attribute the strongest observed line at 2.15 microns to tellurium (Te) III from an M1 line list of heavy elements presented in [45], where the line wavelengths are experimentally calibrated according to the NIST database [156]. Te belongs to the second r-process peak. With the M1 line list, we model kilonova emission line spectra under the assumption that photons from forbidden lines produced outside the photosphere freely escape from the ejecta. The collision strengths of Te III are taken from an R-matrix calculation [43] and those of other ions are obtained by using an atomic structure code HULLAC [157]. The abundance pattern is chosen to be the solar r-process but we separate "light" and "heavy" elements at an atomic mass of 85 and introduce a parameter, the abundance ratio of the two (see figure 14). The ionization fractions are fixed to be (Y +1 , Y +2 , Y +3 ) = (0.2, 0.5, 0.3) motivated by the Te ionization evolution in kilonova ejecta [41]. The line shape is approximated by a Gaussian with a line broadening velocity of 0.08c, which is the same as the photospheric velocity. The mass in the line forming region is estimated by assuming that the observed line luminosity, 5 × 10 38 erg s −1 , is locally generated by radioactivity of r-process nuclei, corresponding to ∼ 0.02M . Given the abundance pattern and ionization state, the mass of Te III in the line forming region is ≈ 8 · 10 −4 M . The electron temperature of the line forming region is then determined such that the total line luminosity agrees with the observed one. The estimated electron temperature is ∼ 3000 K, which is slightly higher than that derived from the pure neodymium nebular modeling [158]. This is because the cooling by tellurium ions is more efficient than neodymium.
We find that [Te III] 2.10 µm line is indeed the most outstanding emission line around 2 microns. Several weaker lines also contribute to the flux around 3-4 microns. There is another potential line feature around 4.5 microns in the NIRSpec spectrum. The location of this feature is consistent with [Se III] 4.55 µm and [W III] 4.43 µm as pointed out by [45] for the kilonova AT2017gfo. From the spectral modeling, we obtain the total ejecta mass of ∼ 0.05-0.1M , which agrees with the one obtained from the light curve modeling ∼ 0.1M .
Here we show a brief estimate of the Te III mass from the observed line at 2.15 microns ( 3 P 0 -3 P 1 ). The collisional excitation rate per Te III ion from the ground level ( 3 P 0 ) to the first exited level ( 3 P 1 ) is given by where n e and T e are the thermal electron density and temperature, Ω 01 ≈ 5.8 is the collision strength [43], E 01 ≈ 0.6 eV is the excitation energy, g 0 is the statistical weight of the ground level. Assuming that the ejecta mass in the line forming region is 0.02M expanding with 0.08c and the ions are typically doubly ionised, we estimate n e ∼ 3 · 10 5 cm −3 , and thus, the line emissivity per Te III ion is 10 ≈ 2.5 · 10 −14 n e 3 · 10 5 cm −3 erg/s, (5) where T e = 3000 K is used. Combining the line emissivity with the observed line luminosity in 2.25 ± 0.23 µm, L line ≈ 3 · 10 38 erg/s, we obtain The mass estimated from the line is somewhat dependent on T e and n e . However, we emphasise that, with T e ≈ 3000 K and n e ≈ 3 · 10 5 cm −3 , the line luminosity is consistent with the radioactive power in the line forming region. It is also interesting to note that the Te III mass of 10 −3 M is in good agreement with the one obtained based on the same line seen in AT2017gfo at 10.5 day [42]. While we conclude that the observed line feature at 2.1 microns is most likely attributed to Te III, it is important to note that there are caveats associated with our modeling. One obvious caveat is that the model does not include E1 lines. Lanthanides and actinides have E1 transitions between low-lying levels in the mid-IR [discussed in 37]. Due to their lower abundances, these lines are expected to be weaker compared to the Te line if collisional excitation dominates the excitation processes. However, as we make an implicit assumption that their E1 lines contribute to the opacity in the mid-IR, they may produce P-Cygni like features, see, e.g., [37,159]. For example, Ce III has a strong line at 2.07 µm with log gf = −1.67 [159]. We estimate that its line optical depth is 0.1 at 700 K with ∼ 0.05M and ∼ 0.1c even if Ce is purely in Ce III. However, more careful analyses including non-LTE effects are needed to quantify it. Another caveat is that the opacity of lanthanides is expected to have some wavelength dependence. Including this effect may also affect the spectral modelings.

Alternative progenitor possibilities
Our interpretation of GRB 230307A provides a self-consistent model for the source in which the temporal and spectral evolution, as well as the source location, can be readily explained. The kilonova has marked similarities with AT2017gfo providing a robust indication of its origin, and we do not need to postulate new and unseen phenomena to explain it. However, it is also relevant to consider alternative possibilities. In particular, given the location of the galaxy at z = 3.87, it is important to consider if the burst could originate at that redshift.

GRB 230307A as a high redshift, highly energetic GRB
The nearby galaxy H1 (F277W(AB)=27.5± 0.1 , r proj = 0.3 arcsec) with a spectroscopic redshift of z=3.87 has a relatively low probability of occurring by chance (∼ 5−10%, see section 4.3). This galaxy has a comparable P chance to G1 (the z = 0.065 galaxy). The host-normalised offset for H1 is ∼ 2.5, which is large but not unprecedented for long GRBs [116]. However, assuming the late time light at the GRB position is all from the transient, it does not lie on the stellar field of this galaxy, which is unusual, for example, in the samples of [116,160,161], there is only one (of > 100) sub-arcsecond localised GRB not on the stellar field of its host.
At z = 3.87, the inferred isotropic energy release and luminosity of GRB 230307A would be E iso = 1.2 × 10 56 erg and L iso = 1.7 × 10 56 erg s −1 (using a 64 ms peak flux). This is approximately an order of magnitude more energetic and two orders of magnitude more luminous than any other previously identified GRB [52].
If at z = 3.87, we can have some confidence that GRB 230307A would be the most energetic burst ever detected by Swift or Fermi, including those without redshift or even afterglow identifications. In Table 4, we tabulate the most fluent GRBs observed by Fermi. Most of these have either redshifts or optical detections, which constrain z < 6 via the detection of the source in the optical band. This leads to a set of measured or maximum E iso values. For events without any redshift information, we can place a conservative upper redshift limit of z = 16. No GRBs detected by Fermi without a redshift can have energy over 10 56 erg unless they lie beyond z ∼ 20. Hence, GRB 230307A is sufficiently rare, if at z = 3.87, that events like it occur less frequently than once per decade across the Universe (i.e. no more than one in the combined lifetimes of Swift and Fermi).
The energetics of the burst at this redshift lie would lie far beyond those of the general GRB population and beyond those suggested as the upper limit for GRB energetics [162]. The only population of core-collapse GRBs whose energetics have been suggested to approach this value are those from first-generation population III stars [163][164][165]. It is not expected that such stars should exist at z ∼ 4. However, while a pop-III origin may alleviate energy concerns, the properties of the GRB and its optical/IR counterpart do not resemble the predictions for pop-III stars. In particular, pop III GRBs are suggested to have particularly long durations given the mass and radii of their progenitors [165], and so require extremely long durations of engine activity to enable jet breakout. However, the ∼ 35 s duration of GRB 230307A and its rapid variability do not readily fit this expectation.
If one ascribes the GRB to a stellar collapse event, considering the afterglow's properties and associated supernovae is also relevant. Firstly, at 28.5 days, the JWST spectral observations are inconsistent synchrotron emission, suggesting that the counterpart must be dominated by another source in the mid-IR and the earlier K-band points. This excess, which in our preferred model is explained by a kilonova, would have to be due to the supernova or shock breakout if at z = 3.87. The K-band (rest frame B-band) light would reach a peak of M B (AB) < −23.5 on a timeframe of <2-days (rest-frame) before decaying at a rate of > 1 mag day −1 for the next four days, or as a power-law decay, a rate of approximately t −4 . This appears too rapid for radioactively powered transients, at least based on the sample to date (note that at z = 3.87, the timescales are a factor of ∼ 5 faster than in the z = 0.065 scenario due to cosmological time dilation). The most likely option for such emission would be shock-breakout, which may begin blue but rapidly cool. There are simulations for the shock breakout associated with pop-III supernovae, which show an early peak [166]. However, this emission peaks in the UV to soft X-ray regime and at luminosities below that seen in GRB 230307A. Indeed, taking 28.5 days as a baseline; for a plausible maximum Pop III radius (e.g. 2000 R , [167]) and a luminosity of L bol ∼ 10 43 erg s −1 the inferred temperature is T ∼ 300, 000 K. This is incompatible with the spectral shape seen in the counterpart to GRB 230307A which peaks at > 4.5 microns (T < 3000K at z = 3.87). Dust or metal line blanketing could alleviate this discrepancy to some degree, but it would be extreme to explain the observations. It would also come at the cost of an even higher intrinsic luminosity. Conversely, the radius at which the luminosity and temperature would be consistent is extremely large (∼ 0.3 pc) and indeed would require super-luminal expansion to reach from a single explosion within the time since burst. These constraints become even more extreme for the Kband observations at 11.5 days, where the luminosity is > 50 times higher. However, we lack detailed information regarding the spectral shape at this time. We can conclude that a thermal transient launched at the time of GRB 230307A cannot explain the observed source at z = 3.87.
We should consider if GRB 230307A could be related to an explosion which bears little to no similarity to long-GRB progenitors. Given the inferred energetics, this is not an unreasonable proposition. However, the emission is too bright and too fast for, for example, the fast blue optical transients (e.g. [168]), the fastest of which have halftimes of ∼ 4 days [169, c.f. < 1 day for GRB 230307A at z = 3.87]. A further alternative may be a relativistic tidal disruption event. This would face significant challenges with the rapid variability timescales seen in the prompt emission and the non-nuclearity of the source within the galaxy at z = 3.87. Putting aside these concerns, the peak optical/IR luminosity is comparable to AT2022cmc [170,171], but the evolution is too rapid and the dynamic range too large.
Finally, it is possible that the red excess seen at later times is not directly related to the progenitor or the transient but is a result of the re-processing of the GRB radiation by material within the host galaxy. In particular, for GRB 211211A [172] suggest that an alternative explanation for the emission could be the heating of dust. However, this model also encounters significant issues at z = 3.87. In particular, the observed K-band excess is a rest-frame B-band excess, much bluer than expected for dust heating. If this represented the peak of the thermal spectrum, it would be above the sublimation temperature of the dust. Alternatively, if it were the blue tail of a much cooler black body, the luminosity would be extremely high.
A final challenge to the high−z scenario is that the afterglow is detected in the UVOT-white and ground-based g-bands. These observations all have substantial sensitivity blue-ward of Lyα at z=3.87. A typical column from the intergalactic medium should attenuate ∼ 50% of the light in these bands, inconsistent with observations. Indeed, for a typical β = 1 spectrum, we expect to observe white − i = 2.8 and g − i = 1.9, approximately 3 and 4 σ away from the observed colours. There is significant variation in the absorption strength as a function of the line of sight, so a low (or near zero) absorption column would alleviate this tension. The sample of [173] implies that at z ∼ 4, perhaps 10% of galaxies have such low absorption sight lines.
Hence, while the proximate galaxy could indicate a high redshift, there are few other indications in the transient properties that would support this interpretation. In particular, neither standard thermal nor non-thermal emission can explain the observed counterpart properties. If the burst does arise from z = 3.87, it requires a new kind of explosion, unlike any seen until now. In practice, such explosions could be extremely rare: the volumetric rate of GRBs with E iso > 10 56 erg is minimal, but postulating them is unnecessary when a robust, physically motivated explanation can be obtained for a lower redshift solution.

Other cosmological scenarios
We should also consider a further option which is that GRB 230307A does not reside at either z = 0.065 or z = 3.87 but is a chance super-position with both galaxies. In this scenario, the actual host is undetected or is one of the other galaxies within the field. The absence of direct redshift measurements makes placing constraints on this scenario challenging. However, we can use the non-detection of the late-time JWST magnitudes to limit the brightness of a supernova component at any redshift.
To quantify the exclusion of "normal" long GRBs at intermediate redshift we utilize model light curves for SN 1998bw from MOSFiT [174] calculated at a range of redshift from 0.05 < z < 4.0 (we take z = 4 as an upper limit for the redshift the GRB based on the observed g-band detections together with the detection of continuum emission to 5300Å [13] and similar faint trace seen to ∼5100Åin our X-shooter spectrum). At each redshift, we compare our observed JWST photometry in each band to the model predictions at that time and report the most constraining limit (e.g. the lowest ratio of F obs /F 98bw ). This is shown in Figure 15. At all plausible redshifts, any supernova must be at least a factor ∼ 3 fainter than SN 1998bw at similar times. For any redshift where the burst energetics fall within the range seen in the bulk GRB population (z < 1.2), any supernova must be a factor > 100 fainter than SN 1998bw. Hence, there is no route in which GRB 230307A can arise from a classical long GRB (E iso < 10 55 erg with an associated broad-lined Ic supernova. The strength of this rejection is predominantly based on the faintness of the source in the bluer bands, whereas at the first epoch, the redder bands are substantially brighter than SN 1998bw at z = 3.87. If the burst lies at an intermediate redshift, dust models may become more appealing. In particular, for a low to moderate redshift, the luminosty and timescales may be a suitable match (e.g. for GRB 211211A [172] find plausible explanations at z ∼ 0.5). However, in this case, the lack of a supernova to extremely deep limits would be surprising, as would the non-detection of the host galaxy.

Galactic objects
In the absence of a robust absorption redshift, it is also necessary to consider if GRB 230307A could arise from a Galactic system.
The very faint magnitudes and extreme red colours observed at late times can effectively rule out X-ray binary outbursts. For example, with an M-dwarf companion (absolute magnitude 9), the distance to the source would be ∼ 100 kpc, and larger for any more massive star, while the late-time colours of the source are not stellar. In practice, given that the source is transient, and we may also expect some contribution from an accretion disc, the overall properties cannot be remedied within the accreting binary framework.
The inferred energetics for Galactic systems are E iso ∼ 5 × 10 43 (d/10 kpc) 2 . This energetic output is within the bounds of giant outbursts from magnetars. For example, the giant flare of SGR 1806-20 had an inferred isotropic energy of E iso ∼ 2 × 10 46 erg (see [175], although subsequent downward revisions in its distance lower this somewhat [176]). However, GRB 230307A does not appear to meet the requirements of a magnetar. In particular, it is at high Galactic latitude (-36 degrees) and far from any plausible star formation that could give rise to a young neutron star. Furthermore, the emission in all magnetar outbursts is dominated by a very short pulse followed by decaying emission in which the pulse period of the neutron star is visible. This is not the case for GRB 230307A.
GRBs have previously been suggested to arise from the tidal disruption and accretion of rocky material onto a neutron star [177], and such events are seen in the case of white dwarfs [178]. However, for accretion onto a neutron star, we would usually expect to observe a relatively soft outburst (e.g. in the model of [177] the temperature is ∼ 10 keV). The spectrum we observe for GRB 230307A consists of evolving synchrotron emission which does not contain a thermal component and is inconsistent with directly observing heating accreting material. Indeed, even for near direct accretion, it is not clear how such a spectrum would be formed in an accreting neutron star scenario. Indeed, in this case, the evolution to very low temperatures on a timescale of ∼ 60 days would also not be natural.
Hence, we conclude that no known Galactic systems could explain the observed properties of GRB 230307A.

GRB 230307A as a white-dwarf -neutron star merger
A final alternative is that GRB 230307A is related to the merger of a white dwarf and a neutron star. Although this is still a "compact object merger", such mergers are very different from those of neutron stars with another neutron star or a black hole. In particular, simulations show that no r-process material is produced [179], and so we should not expect the very red emission. Although there are suggestions that GRB 211211A could have been produced by a WD-NS merger [19,59], it is unclear if these could readily explain the detailed spectrophotometric evolution of GRB 230307A.
White dwarf neutron star mergers are appealing, because the long-duration of the gamma-ray emission could suggest a less compact remnant merger event. The wider separation of the binary at the disruption of the white dwarf produces disk accretions times from 100 to 1000s, matching the long duration of these bursts [180]. However, we expect the accretion rate in the mergers to be low, producing less-powerful, and hence, less luminous GRBs. Current WD/NS merger simulations predict a range of light curves that span the emission from GRB 230307A. A Ca feature does exist that is close to the observed line feature, but the models are too blue to explain the shape of the spectra [179]. This is because WD/NS mergers do not produce elements much heavier than the iron peak elements. As we have noticed in matching kilonova models, these elements do not have strong lines beyond ∼ 20,000Å. The subsequent emission above these wavelengths is very weak in these models. In addition, WD/NS mergers are expected to have fairly week kicks to ensure that the binary remains bound and these mergers are expected to have much lower offsets than neutron star mergers. However, the mass of the best-candidate host galaxy is sufficiently low that the observed offset for this burst can be attained [181].    Table 4 The properties of the brightest 10 bursts observed by Fermi/GBM. Eight of these have afterglow identifications which place limits on the redshift, while it is likely that all originate from z < 16. As expected for the (observationally) brightest bursts there is a preference for low redshift, with the bursts all arising at much less than the mean for Swift GRBs [182]. In 3-4 cases a supernova has been identified, meaning the bursts arise from core collapse (with GRB 221009A still ambiguous at the time of writing). In 2 cases a kilonova origin appears more likely, while for the remaining 4 the absence of afterglows or the redshifts render such diagnostics impossible. However, it would appear that samples of bright, long-GRBs may contain significant number of compact object mergers.   Table 6 Light curve model parameters, priors and posterior medians and 90% credible intervals.

Author contribution statements
AJL led the project, including the location of the afterglow and kilonova and the JWST observations. BPG first identified the source as a likely compact object merger, was co-PI of the Chandra observations, and contributed to analysis and writing. OS contributed to afterglow and kilonova modelling and led the writing of these sections. MB was involved in kilonova modelling, EB contributed to interpretation, placing the burst in context and high energy properties. KH was involved in kilonova spectral modelling and identified the 2.15 micron feature. LI reduced VLT/MUSE and X-shooter observations and led the host analysis. GPL contributed to afterglow and kilonova modelling. ARE analysed the Chandra observations, BS reduced and analysed VLT observations. NS contributed to afterglow and kilonova modelling. SS was responsible for placing the burst afterglow in context and demonstrating its faintness. NRT contributed to analysis, interpretation and writing. KA was involved in the ULTRA-CAM observations and interpretation. GA led the ATCA observations. GB reduced the JWST NIRCAM data. LC processed and analysed the MUSE observations. VSD is the ULTRACAM PI. JPUF contributed to the interpretation. WF was the PI on the Chandra observations and contributed to discussion. CF contributed to the theoretical interpretation. NG was involved in host analysis. KEH, GP, AR, SDV, SC, PDA, DH, MDP, CCT, AdUP and DA contributed to ESO observations and discussion. DW contributed to spectral and progenitor modelling and discussion. MJD, PK, SP, JM, SGP, IP, DIS contributed to the ULTRACAM observations. GL investigated potential similarities with other transients. AT, PAE, BS and JAK contributed to the Swift observations. MF extracted and flux-calibrated the TESS light curve. SJS analysed the JWST spectral lines, and contributed to interpretation and writing. HFS performed the BPASS-hoki-ppxf fits to the integrated MUSE flux and contributed the associate figure and text. All authors contributed to manuscript preparation through contributions to concept development, discussion and text.

Acknowledgements
We dedicate this paper to David Alexander Kann, who passed on March 10. The final messages he sent were regarding follow-up of GRB 230307A, and we hope it would satisfy his curiosity to know the final conclusions. This work is based on observations made with the NASA/ESA/CSA James Webb Space Telescope. The data were obtained from the Mikulski Archive for Space Telescopes at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-03127 for JWST. These observations are associated with program #4434 and 4445. Support for Program numbers 4434 and 4445 was provided through grants from the STScI under NASA contract NAS5-03127.
This paper is partly based on observations collected at the European Southern Observatory under ESO programme 110.24CF (PI Tanvir)

Data availability
JWST data are directly available from the MAST archive. Chandra and Swift data are also in the public domain. ESO and Gemini data are stored in their respective archives and will be available to all once the proprietary period expires. Data can be obtained from the corresponding author between the date of publication and the end of the proprietary period. This research has made use of Fermi data which are publicly available and can be obtained through the High Energy Astrophysics Science Archive Research Center (HEASARC) website at https://heasarc.gsfc.nasa.gov/W3Browse/ fermi/fermigbrst.html

Code availability
Much analysis for this paper has been undertaken with publically available codes and the details required to reproduce the analysis are contained within the manuscript. Fig. 6 The distribution of measure fluence from the Fermi/GBM catalogue [105]. The solid line shows an expected slope of −3/2 for a uniform distribution. The faint end deviates from this line because of incompleteness. At the brighter end, there are three bursts which appear to be extremely rare, GRB 130427A, GRB 221009A [52] and GRB 230307A. To indicate the apparent rarity we also plot lines representing the expected frequency of events under the assumption of a −3/2 slope. We would expect to observe bursts akin to GRB 230307A only once per several decades. Fig. 7 A comparison between the prompt fluence (in the 15-150 keV band) and the X-ray flux at 11 hours for Swift GRBs with GRB 221009A and GRB 230307A added, updated from [106,107]. The general 1:1 trend between the prompt fluence and X-ray brightness can clearly be seen, although it has a significant scatter, although a very rare event, GRB 221009A apparently lies on the same relation. However, GRB 211211A and GRB 230307A are clearly outliers to this relation with GRB 230307A occupying a region devoid of other GRBs. This very faint afterglow compared to the prompt emission may be related to a location at large projected offset from its host galaxy in which the density of the ambient medium is very low.            15 Limits on supernova similar to SN 1998bw as a function of redshift, based on the most constraining detection with JWST. At any redshift for which GRB 230307A would not be the most energetic GRB ever observed, any supernova is at least a factor of ∼ 100 fainter than SN 1998bw.