We present a summary of gravitational-wave (GW) follow-up using the Las Cumbres Observatory global network of telescopes during the third (O3) and fourth (O4) observing runs of the GW detectors. As in O2, we implemented the Gehrels et al. 2016 galaxy-targeted strategy. Here we test its efficacy in O3 and O4 and analyze the Las Cumbres Observatory response time and depth for nine GW alerts that showed a possibility of having an electromagnetic counterpart (GW190425, GW190426_152155, S190510g, GW190728_064510, GW190814, S190822c, GW191216_213338, S240422ed and S250206dm). We find that Las Cumbres Observatory is able to begin observations in response to GW alerts within minutes of the alert, with the observations being deep enough to detect possible GW170817-like kilonovae out to a median distance of 250 Mpc. In this sense a global rapid-response network of telescopes like Las Cumbres is an excellent GW follow-up facility. However, the galaxy-targeted follow-up strategy was much less efficient in O3 and O4 than originally predicted, given the larger than assumed GW localizations. We conclude that coordination between various facilities to include both wide-field and rapid-response capabilities is required to achieve efficient and comprehensive follow-up of GW events.

Quasi-periodic eruptions (QPEs) are recurring bursts of X-ray radiation originating from supermassive black holes (SMBHs). They are an unprecedented type of structured, high-amplitude SMBH variability, but the physical origins of their regularity, timescales, energetics, and emission are uncertain. We present new XMM-Newton observations of the QPEs in ZTF19acnskyy/“Ansky”, constituting the deepest observations of individual bursts in any source thus far. The X-ray spectra reveal time-evolving P Cygni profiles comprising blueshifted absorption and redshifted emission from L-shell transitions of Fe XIX-XXIV, with column densities N_H∼ 10^22-23 cm^-2 and bulk velocities of |v_w/c|∼ 0.2, indicating relativistic mass ejections during each eruption. We construct a time-dependent analytical model of a wind turning on to self-consistently compute its evolving luminosity and ionization properties, and find that the light curve and spectral lines can be simultaneously produced by a wide-angle outflow with Ṁ∼ 10^-9-10^-8 M_⊙ s^-1 kinetically powering the X-rays with an efficiency of L_X/Ė_K∼ 0.1. Each eruption ejects ∼ 10^-3 M_⊙ and ≳ 10^49 erg of kinetic energy, setting an upper bound on the QPE lifetime of ≲30 years if the underlying mass reservoir is ∼1 M_⊙, and implying that the bursts may result in detectable multiwavelength signatures of reverberation and feedback. These measurements provide new quantitative constraints on QPE energetics, emission mechanisms, and the mass/energy they recycle into their circumnuclear environments, as well as an observational probe for direct comparison with physical models and hydrodynamical simulations of QPEs.

Optical tidal disruption events (TDEs) exhibit extremely broad emission lines (≈ 10^3-10^4 km s^-1) and are observationally classified into four spectroscopic types: H-dominated, He-dominated, H+He, and featureless. The prevalent H+He class often displays Bowen fluorescence lines (notably and ), features that are rarely observed in active galactic nuclei and whose origin has remained poorly understood. We present the first unified radiative transfer framework that reproduces all four TDE spectroscopic classes using simulations of optically thick, outflowing envelopes with solar composition. Our models successfully capture both the continuum properties and key spectral features, including strong , and Bowen emissions. We demonstrate that the spectroscopic diversity of TDEs is primarily governed by the gas ionization state, controlled by the ratio of injected luminosity to envelope mass. As the ionization level decreases, the observed sequence of spectroscopic classes emerges naturally, transitioning from featureless to He-dominated, to Bowen-dominated, and finally to H-dominated spectra. We further show that electron scattering in the optically thick outflow is the dominant mechanism responsible for the extreme line widths, linking line profiles directly to the physical properties of the wind. The model also explains the observed correlations with luminosity, black hole mass, and the relative stability of spectral classifications during TDE evolution. This work establishes a unified physical framework for TDE spectroscopy, providing new insight into the emission mechanisms, energetics, and outflow structure of these transient events, and offering a practical pathway for interpreting and fitting observed spectra.

We predict the gamma-ray line emission from r-process nuclei synthesized in the ejecta of the accretion-induced collapse (AIC) of a magnetized, rapidly rotating white dwarf. Using ejecta from a two-dimensional general-relativistic neutrino-magnetohydrodynamic simulation, further evolved with a radiation-hydrodynamics code coupled to an in situ nuclear reaction network, we construct angle-dependent gamma-ray spectra in the 0.01-10 MeV band via composition-dependent ray tracing through the ejecta. The emission between 1 and 10 d is dominated by 132I (t1/2 = 2.3 h), continuously replenished by the decay of its parent 132Te (t1/2 = 3.2 d), with additional contributions from 131I, 133Xe, and 132Te. At t greater than or similar to 20 d, 56Co (from 56Ni decay) becomes the primary emitter. The simultaneous presence of r process and iron-peak gamma-ray lines is distinctive of AIC ejecta and absent in binary neutron star mergers, where iron-peak nuclei are generally not synthesized. Comparing with the 3 sigma continuum sensitivities of planned MeV gamma-ray telescopes (COSI, AMEGO-X, e-ASTROGAM, GRAMS, GammaTPC), we find the brightest r-process lines detectable to 10 Mpc by GammaTPC and GRAMS, with the signal approaching their sensitivity threshold at 30 Mpc. The r-process spectral features survive time integration over 30 d exposures, demonstrating robustness against the long observation times required by gamma-ray detectors.

The source of the optical/UV emission in tidal disruption events (TDEs) remains an enduring question in the field. Connecting the observed emission to the source is critical for both our understanding of these transients and for using TDEs to study the efficiency of super-Eddington accretion and black hole growth. To explore this connection, we ran time-dependent 1D radiation hydrodynamic simulations of TDE emission with the Sedona Monte Carlo radiative transfer code, focusing on the reprocessing paradigm. Our simulations follow a compact, evolving X-ray and EUV bright source and surrounding reprocessing outflow over multiple months, using luminosities and mass flow rates consistent with hydro simulations of tidal disruptions. We determine the efficiency of reprocessing as a function of time in this dramatically changing environment and reproduce key observables, including timescales, luminosities, and color evolution. Notably, we see a strong wavelength dependence in the emission timescale due to reprocessing effects. Early on, there is an X-ray flare that quickly fades as material builds up and obscures the hot source. At the same time, the optical/UV luminosity begins to rise. Though the optical/UV light curve has a similar shape to the bolometric light curve, the optical peak is offset by similar to 3 weeks from the bolometric peak due to the time required to build up the reprocessing layer. This implies that early time, high-energy emission may be missed for TDEs discovered in optical surveys, and the initial disruption and mass return time to the black hole may occur earlier than optical light curves suggest.

Featureless optical and ultraviolet (UV) spectra are a puzzling signature to emerge from recent observations of luminous fast blue optical transients (LFBOTs) and some tidal disruption events (TDEs). We describe the landscape of source and gas properties that are expected to form H, He I and He II emission lines, and map spectral types to the parameter space of luminosity and system radius. Using one-dimensional radiative transfer calculations, we show that high source luminosities (L > 10^44 erg s^-1) and compact ejecta radii (r < 10^14 cm) produce featureless spectra due to the high temperature and ionization state of the emitting medium. Intermediate luminosities and moderately compact systems can generate He II-dominated spectra, while lower luminosities and more extended atmospheres result in conspicuous H and He I emission. Large expansion velocities (v ≥ 0.1c) can further broaden lines such that they blend into the continuum. Featureless UV spectra may require even more extreme ionization conditions or velocities to suppress the many intrinsically strong metal lines at those wavelengths. Applying this framework to understand the absence of features observed in LFBOTs and featureless TDEs, we find that non-homologous, compact outflows are likely necessary for featurelessness to persist in optical and UV spectra.

Neutron star mergers are a leading site of r-process, producing radioactively powered optical and infrared transients known as kilonovae. Observations of the kilonovae AT2017gfo, associated with the gravitational-wave event GW170817, and AT2023vfi, associated with GRB 230307A, have enabled measurements of the mass of ejected r-process material and the identification of heavy elements in the ejecta. However, late-time observations reveal strong infrared emission with temperature below 1000 K, which is difficult to explain by atomic absorption and emission processes alone. In this paper, we show that kilonova ejecta provide conditions favorable for the formation of dust grains composed of refractory r-process elements including Zr, W, and Os. We calculate the kinetic formation of dust grains using reaction rate coefficients of W as a proxy, finding that dust forms efficiently, particularly in slow ejecta. This stands in contrast to a previous study that relied on a classical nucleation framework. By performing radiative transfer simulations that incorporate dust formation, we demonstrate that r-process dust naturally explains the observed late-time infrared emission. The formation and abundance of r-process dust are highly sensitive to the ejecta mass, composition, and expansion velocity. Infrared emission from r-process dust can therefore serve a new probe of heavy-element production in neutron star mergers.

A long-lived central engine embedded in expanding supernova ejecta can alter the dynamics and observational signatures of the event, producing an unusually luminous, energetic, and/or rapidly-evolving transient. We use two-dimensional hydrodynamics simulations to study the effect of a central energy source, varying the amount, rate, and isotropy of the energy deposition. We post-process the results with a time-dependent Monte Carlo radiation transport code to extract observational signatures. The engine excavates a bubble at the centre of the ejecta, which becomes Rayleigh-Taylor unstable. Sufficiently powerful engines are able to break through the edge of the bubble and accelerate, shred, and compositionally mix the entire ejecta. The breakout of the engine-driven wind occurs at distinct rupture points, and the outflowing high-velocity gas may eventually give rise to radio emission. The dynamical impact of the engine leads to faster rising optical light curves, with photon escape facilitated by the faster expansion of the ejecta and the opening of low-density channels. For models with strong engines, the spectra are initially hot and featureless, but later evolve to resemble those of broad-line Ic supernovae. Under certain conditions, line emission from ionized, low-velocity material near the centre of the ejecta may be able to escape and produce narrow emission similar to that seen in interacting supernovae. We discuss how variability in the engine energy reservoir and injection rate could give rise to a heterogeneous set of events spanning multiple observational classes, including the fast blue optical transients, broad-line Ic supernovae, and superluminous supernovae.

Neutron star-bearing compact-object mergers can create a kilonova, a transient powered by the radioactive decay of material synthesized via rapid neutron capture (r-process). During the merger, a fraction of the ejecta is launched at high speeds (≳ 0.4c) that can lead to neutrons evading capture onto seed nuclei, resulting in free neutrons. The free neutrons then decay and inject additional energy that power early-time (≲1 day) emission. Here, using , we present a grid of the first non-local thermodynamic equilibrium and frequency-dependent opacity radiative transfer simulations to predict early X-ray/UV/optical emission from free neutron decay that span the expected theoretical range of free neutron mass, mixing with r-process material, and extent of the high-velocity ejecta tail. The emission can be characterized by an SED that rapidly shifts from X-rays of ∼few×10^41-10^42 erg s^-1 in the first ∼ minutes to a far-UV and near-UV peak on the timescale of ∼ minutes to hours, followed by enhanced optical and IR emission for sufficiently large free neutron masses. The properties of the free neutron ejecta are most distinguishable, in principle, at extreme-UV and far-UV wavelengths, with SED peak flux and wavelength determined by the maximum velocity of the ejecta and the mixing. In the far-UV and near-UV, free neutron ejecta masses as small as 10^-7 M_⊙ exhibit a unique bump compared to neutron-free models, demonstrating the importance of upcoming UV missions like UVEX and ULTRASAT to constraining the nucleosynthetic environment of neutron star mergers.

We present SEDONA-GesaRaT, a rapid code for supernova radiative transfer simulation developed based on the Monte Carlo radiative transfer code SEDONA. We use an atomic physics neural network, which has a convolutional neural network structure, to solve for the nonlocal thermodynamic equilibrium (NLTE) atomic physics level population calculation. SEDONA-GesaRaT is trained and validated on 119 1D type Ia supernova (SN Ia) radiative transfer simulation results. It has been applied to the 3D SN Ia explosion model N100 to perform a 3D NLTE radiative transfer calculation. The spatially resolved linear polarization data cubes of the N100 model are successfully retrieved with a high signal-to-noise ratio using the integral-based technique. The 3D LTE and NLTE radiative transfer calculations have been validated by comparisons to SEDONA program simulation results. The typical computation cost of a new 3D NLTE spectropolarimetry simulation using SEDONA-GesaRaT is only ∼3000 core-hours, which will be a factor of ∼100 faster relative to a SEDONA simulation. The large computational speed-up factor makes SEDONA-GesaRaT promising for future large-scale simulations that systematically study the internal structures of SNe. However, the ultimate advantages of SEDONA-GesaRaT will depend on the required computational cost of retraining to accommodate improvements in physical realism in radiative transfer and a wider range of SN models.

We present an extensive photometric and spectroscopic ultraviolet–optical–infrared campaign on the luminous fast blue optical transient (LFBOT) AT 2024wpp over the first ∼100 days. AT 2024wpp is the most luminous LFBOT discovered to date, with L _pk ≈ (2–4) × 10 ^45 erg s ^−1 (5–10 times that of the prototypical AT 2018cow). This extreme luminosity enabled the acquisition of the most detailed LFBOT UV light curve thus far. In the first ∼45 days, AT 2024wpp radiated >10 ^51 erg, surpassing AT 2018cow by an order of magnitude and requiring a power source beyond the radioactive ^56 Ni decay of traditional supernovae. Like AT 2018cow, the UV–optical spectrum of AT 2024wpp is dominated by a persistently blue thermal continuum throughout our monitoring, with blackbody parameters at a peak of T > 30,000 K and R _BB / t ≈ 0.2 c –0.3 c . We find evidence for cooling until ∼10 days; thereafter, T ≳ 20,000 K is maintained. We interpret the featureless spectra as a consequence of continuous energy injection from a central source of high-energy emission that maintains high ejecta ionization. After 35 days, faint (equivalent width (EW) ≲ 10 Å) H and He spectral features with kinematically separate velocity components centered at 0 and −6400 km s ^−1 emerge, implying spherical symmetry deviations. A near-infrared excess of emission above the optical blackbody emerges between 20 and 30 days, with a power-law spectrum F _ν _,NIR ∝ ν ^−0.3 at 30 days. We interpret this distinct emission component as either reprocessing of early UV emission in a dust echo or free–free emission in an extended medium above the optical photosphere. LFBOT asphericity and multiple outflow components (including mildly relativistic ejecta), together with the large radiated energy, are naturally realized by super-Eddington accretion disks around neutron stars or black holes and their outflows.

We present the first end-to-end calculation connecting the accretion-induced collapse (AIC) of a magnetized, rapidly rotating white dwarf to observable kilonova signatures, combining two-dimensional (2D) general-relativistic neutrinomagnetohydrodynamic simulations, followed by radiation hydrodynamics with in-situ nuclear network and 2D Monte Carlo radiative transfer with spatially resolved heating rates. Unlike all previous unmagnetized AIC models-which predicted proton-rich, Ni-56-dominated ejecta-strong magnetic fields eject ti 0 . 2 M-circle star of neutron-rich material ( ( Y-e ) similar to 0.24) on dynamical time-scales, before neutrino irradiation can raise the electron fraction, enabling strong r-process nucleosynthesis up to and beyond the third peak. The resulting kilonova is lanthanide-rich (X-lan approximate to 8 per cent ) and dominated by near-infrared emission. We compute synthetic light curves in the Large Synoptic Survey Telescope and James Webb Space Telescope bands and find striking agreement, without parameter tuning, between the observations of AT 2023vfi/GRB 230307A and our broadband light curves for polar viewing angles. These results establish magnetized AIC as a viable channel for heavy r-process element production and a compelling progenitor candidate for long-duration gamma-ray bursts with kilonova signatures.

We present the first three-dimensional study of the asymptotic ejecta distributions for a suite of theoretical Type IIp supernovae originating from red supergiant (RSG) progenitors. We simulate using the radiation hydrodynamic code FORNAX from the core bounce through to the first seconds of the neutrino-driven explosion and then follow using a hydrodynamic variant of the code FLASH until the shock breakout of the star and through to homologous expansion of the ejecta into the circumstellar environment. Our studied progenitor models range from 9 to 25 M circle dot, with explosion energies spanning similar to 0.1-1 Bethe. The shock breakout times span the range of similar to 1-4 days, with a breakout time spread by direction ranging from hours to over a day. We find that the dipole orientation of the 56Ni ejecta is well preserved from the first seconds out to the shock breakout. The 56Ni ejecta penetrates through the initially outer oxygen shell, and its global structure is imprinted with small-scale clumping as the ejecta evolve through the stellar envelope. For the majority of our models, the neutron star kick is anti-aligned with the 56Ni ejecta. Models with strongly dipolar ejecta morphology and a massive hydrogen/helium envelope with an inner boundary located deep see as much as similar to 70% of the 56Ni ejecta mixed into that outer envelope, reaching asymptotic velocities ranging from similar to 350 to 3200 km s-1. Supernovae arising from RSG progenitors and exhibiting prominent nickel features generally display significant 56Ni mixing into the stellar envelope.

In order to better connect core-collapse supernova (CCSN) theory with its observational signatures, we have developed a simulation pipeline from the onset of the core collapse to beyond shock breakout from the stellar envelope. Using this framework, we present a 3D simulation study from 5 s to over 5 days following the evolution of a 17 M _⊙ progenitor, exploding with ∼10 ^51 erg of energy and ∼0.1 M _⊙ of ^56 Ni ejecta. The early explosion is highly asymmetric, expanding most prominently along the southern hemisphere. This early asymmetry is preserved to shock breakout, ∼1 day later. Breakout itself evinces strong angle-dependence, with as much as 1 day delay in the shock breakout by direction. The nickel ejecta closely tail the forward shock, with velocities at the breakout as high as ∼7000 km s ^−1 . A delayed reverse shock forming at the H/He interface on hour timescales leads to the formation of Rayleigh–Taylor instabilities, fast-moving nickel bullets, and almost complete mixing of the metal core into the hydrogen envelope. For the first time, we illustrate the angle-dependent emergent broadband and bolometric light curves from simulations evolved in 3D in entirety, continuing through hydrodynamic shock breakout from a CCSN model of a massive stellar progenitor evolved with detailed, late-time neutrino microphysics and transport. Our case study of a single progenitor underscores that 3D simulations generically produce the cornucopia of observed asymmetries and features in CCSNe observations, while establishing the methodology to study this problem in breadth.

We investigate the kilonova emission resulting from outflows produced in a 3D general-relativistic magnetohydrodynamic (GRMHD) simulation of a hypermassive neutron star (HMNS) remnant. We map the outflows into the flash hydrodynamics code to model their expansion in axisymmetry, and study the effects of employing different r-process heating rates. Except for the highest heating rate prescription, we find no significant differences with respect to overall ejecta dynamics and morphology compared to the simulation without heating. Once homologous expansion is attained, typically after similar to 2 s for these ejecta, we map the outflows to the sedona radiative transfer code and compute the spectral evolution of the kilonova and broad-band light curves in various Legacy Survey of Space and Time (LSST) bands. The kilonova properties depend on the remnant lifetime, with peak luminosities and peak time-scales increasing for longer lived remnants that produce more massive ejecta. For all models, there is a strong dependence of both the bolometric and broad-band light curves on the viewing angle. While the short-lived (12 ms) remnant produces higher luminosities when viewed from angles closer to the pole, longer lived remnants (240 ms and 2.5 s) are more luminous when viewed from angles closer to the equator. Our results highlight the importance of self-consistent, long-term modelling of merger ejecta, and taking viewing-angle dependence into account when interpreting observed kilonova light curves. We find that magnetized outflows from an HMNS - if it survives long enough - could explain blue kilonovae, such as the blue emission seen in AT2017gfo.

We present an integral-based technique (IBT) algorithm to accelerate supernova (SN) radiative transfer calculations. The algorithm utilizes “integral packets”, which are calculated by the path integral of the Monte-Carlo energy packets, to synthesize the observed spectropolarimetric signal at a given viewing direction in a 3-D time-dependent radiative transfer program. Compared to the event-based technique (EBT) proposed by (Bulla et al. 2015), our algorithm significantly reduces the computation time and increases the Monte-Carlo signal-to-noise ratio. Using a 1-D spherical symmetric type Ia supernova (SN Ia) ejecta model DDC10 and its derived 3-D model, the IBT algorithm has successfully passed the verification of: (1) spherical symmetry; (2) mirror symmetry; (3) cross comparison on a 3-D SN model with direct-counting technique (DCT) and EBT. Notably, with our algorithm implemented in the 3-D Monte-Carlo radiative transfer code SEDONA, the computation time is faster than EBT by a factor of 10-30, and the signal-to-noise (S/N) ratio is better by a factor of 5-10, with the same number of Monte-Carlo quanta.

We present SEDONA-GesaRaT, a rapid code for supernova radiative transfer simulation developed based on the Monte-Carlo radiative transfer code SEDONA. We use a set of atomic physics neural networks (APNN), an artificial intelligence (AI) solver for the non-local thermodynamic equilibrium (NLTE) atomic physics level population calculation, which is trained and validated on 119 1-D type Ia supernova (SN Ia) radiative transfer simulation results showing great computation speed and accuracy. SEDONA-GesaRaT has been applied to the 3-D SN Ia explosion model N100 to perform a 3-D NLTE radiative transfer calculation. The spatially resolved linear polarization data cubes of the N100 model are successfully retrieved with a high signal-to-noise ratio using the integral-based technique (IBT). The overall computation cost of a 3-D NLTE spectropolarimetry simulation using SEDONA-GesaRaT is only ∼3000 core-hours, while the previous codes could only finish 1-D NLTE simulation, or 3-D local thermodynamic equilibrium (LTE) simulation, with similar computation resources. The excellent computing efficiency allows SEDONA-GesaRaT for future large-scale simulations that systematically study the internal structures of supernovae.

We present X-ray (0.3-79 keV) and radio (0.25-203 GHz) observations of the most luminous fast blue optical transient (LFBOT) AT 2024wpp at z = 0.0868, spanning 2-280 days after first light. AT 2024wpp shows luminous (LX approximate to 1.5 x 1043 erg s-1), variable X-ray emission with a Compton hump peaking at delta t approximate to 50 days. The X-ray spectrum evolves from a soft (F nu proportional to nu-0.6) to an extremely hard state (F nu proportional to nu 1.26) accompanied by a rebrightening at delta t approximate to 50 days. The X-ray emission properties favor an embedded high-energy source shining through asymmetric expanding ejecta. We detect radio emission peaking at L9 GHz approximate to 1.7 x 1029 erg s-1 Hz-1 at delta t approximate to 73 days. The spectral evolution is unprecedented: the early millimeter fluxes rise nearly an order of magnitude during delta t approximate to 17-32 days, followed by a decline in spectral peak fluxes. We model the radio emission as synchrotron radiation from an expanding blast wave interacting with a dense environment ( M similar to 10-3M circle dot yr-1 for vw = 1000 km s-1). The inferred outflow velocities increase from Gamma beta c approximate to 0.07c to 0.42c during delta t approximate to 32-73 days, indicating an accelerating blast wave. We interpret these observations as a shock propagating through a dense shell of radius approximate to 1016 cm and then accelerating into a steep density profile rho CSM(r) proportional to r-3.1. All radio-bright LFBOTs exhibit similar circumstellar medium (CSM) density profiles (rho CSM proportional to r-3), suggesting similar progenitor processes. The X-ray and radio properties favor a progenitor involving super-Eddington accretion onto a compact object launching mildly relativistic disk wind outflows.

Massive stars can end their lives with a successful supernova explosion (leaving behind a neutron star or, more rarely, a black hole), or a failed explosion that leaves behind a black hole. The density structure of the pre-collapse progenitor star already encodes much of the information regarding the outcome and properties of the explosion. However, the complexity of the collapse and subsequent shock expansion phases prevents drawing a straightforward connection between the pre-collapse and post-explosion properties. In order to derive such a connection several explodability studies have been performed in recent years. However, different studies can predict different explosion outcomes. In this article, we show how compactness, which is related to the average density of the star's core, has an important role in determining the efficiency of neutrino heating, and therefore the outcome of the explosion. Commonly, high-compactness progenitors are assumed to yield failed explosions, due to their large mass accretion rates, preventing the shock from expanding. We show by analyzing ∼ 150 2D FLASH and Fornax simulations and 20 3D Fornax simulations that this is not the case. Instead, due to the rapid increase of neutrino heating with compactness, high-compactness progenitors lead to successful shock revival. We also show that 1D+ simulations that include ν-driven convection using a mixing-length theory approach correctly reproduce this trend. Finally, we compare 1D+ models, which we show can reproduce some aspects of multi-D simulations with reasonable accuracy, with other widely used 1D models in the literature.

On 2025 August 18, the LIGO–Virgo–KAGRA collaboration reported gravitational waves from a subthreshold binary neutron star merger. If astrophysical, this event would have a surprisingly low chirp mass, suggesting that at least one neutron star was below a solar mass. The Zwicky Transient Facility mapped the coarse localization and discovered a transient, ZTF 25abjmnps (AT2025ulz), which was spatially and temporally coincident with the gravitational-wave trigger. The first week of follow-up suggested properties reminiscent of a GW170817-like kilonova. Subsequent follow-up suggests properties most similar to a young, stripped-envelope, Type IIb supernova. Although we cannot statistically rule out chance coincidence, we undertake due diligence analysis to explore the possible association between ZTF 25abjmnps and S250818k. Theoretical models have been proposed wherein subsolar neutron star(s) may form (and subsequently merge) via accretion-disk fragmentation or core fission inside a core-collapse supernova—i.e., a “superkilonova.” Here, we qualitatively discuss our multiwavelength dataset in the context of the superkilonova picture. Future higher-significance gravitational-wave detections of subsolar neutron star mergers with extensive electromagnetic follow-up would conclusively resolve this tantalizing multimessenger association.