NOCTURNE. VI. Hour-timescale radio variability in early-stage active galactic nucleus SDSS J122844.81+501751.2

Introduction

Approximately 10% of active galactic nuclei (AGN) can launch and sustain powerful relativistic jets [62], which are often highly variable yet persistently detectable at radio frequencies. Due to their strong radio emission, jetted AGN have often been identified using the radio loudness parameter, $R$¹. Whereas it is still able to pick out the radio-brightest sources, this parameter fails to represent the true diversity of the jet phenomenon in AGN, spanning from powerful relativistic jets to non-collimated outflows and everything in between [62, 43, 63, 33]. Moreover, the $R$ parameter often fails to identify even jetted AGN if they are misaligned or low-power. Alternatively, non-jetted AGN whose hosts exhibit strong star formation can be classified as radio-loud [13, 33, 83].

Radio-loudness faced another hurdle with the discovery of a small sample of extremely variable radio AGN [44]. Most of these sources are narrow-line Seyfert 1 (NLS1) galaxies [16]. The NLS1 identification is based on their optical spectrum: the full-width at half maximum (FWHM) of their broad H$\beta$ emission line is $<2000$ km s$^{-1}$, and their [O III] emission is weak compared to the broad H$\beta$: $S$([O III])/$S$(H$\beta$)$<3$ [60, 22]. A large fraction of them also exhibit strong Fe II emission, confirming the unobstructed view of the central engine. The narrowness of the broad lines can be attributed to the low rotational velocity of the broad-line region (BLR) gas around an intermediate-to-low-mass supermassive black hole ($10^{6}$–$10^{8}M_{\odot}$, [65, 42]). This interpretation is supported by several lines of evidence, such as reverberation mapping studies ([85] (Du et al. 2018)), predominantly turbulence-dominated Lorentzian emission line profiles ([40, 21, 76, 50, 15, 6]), the detection of tidal disruption events in NLS1s ([10, 36, 19]), and the prevalence of disk-like host galaxies with pseudobulges (e.g., [59, 34, 58, 82, 81]). On average, NLS1s accrete at very high rates (average $L/L_{\mathrm{Edd}} \sim0.3$), often exceeding the Eddington limit ([11]), as also indicated by their strong Fe II emission in comparison to the broad H$\beta$ ($R4570$ or $R_{\mathrm{Fe\,II}}$), known to correlate with the Eddington ratio ([49]). These properties suggest that NLS1s are fast-growing AGN in an early evolutionary stage ([51]).

The sample of NLS1s showing extreme variability properties was discovered as part of the Metsähovi Radio Observatory (MRO, Finland) AGN monitoring programme at 37 GHz ([43]), and later all of them have also been detected at 15 GHz by the Owens Valley Radio Observatory (OVRO, USA) ([9, 35]). These sources are intriguing because they are often undetectable, but during flares they can reach Jy-level flux densities. NLS1s were found to host relativistic jets in the late 2000’s ([41, 89, 1]), so the initial hypothesis was that these sources host such jets. One of them was even identified as a gamma-ray emitter and was observed to brighten in X-rays soon after an MRO-detected flare ([70]), supporting the jet hypothesis. However, subsequent medium- and high-resolution radio observations using the Karl G. Jansky Very Large Array (JVLA) and the Very Long Baseline Array (VLBA) did not find any hints of jets ([9, 32, 35]). In VLBA observations at 15 GHz (rms$\sim60$ µJy), none of the sources was detected, ruling out high brightness temperature radio cores seen in other jetted NLS1s, and blazars in general. Multifrequency radio data, including JVLA data from $L$- to $Q$-band, complemented by archival observations down to 144 MHz, revealed that all sources exhibit steep broadband spectra, with several sources undetected at higher frequencies (see Figure 1 in [35]). Only one source, J122844.81+501751.2 (J1228+5017 hereafter), also the focus of this study, showed a flattening of the spectrum at high frequencies, possibly indicating a flat core or a coronal component. The broadband radio spectrum of J1228+5017 is shown in Figure 1, where the slope appears to flatten above 10 GHz. Despite the absence of detectable jets, [35] report blazar-like variability brightness temperatures ($T_{\mathrm{B}} \approx10^{13}$–$10^{15}$ K) and variability Doppler factors ($\delta\approx10$–30) for one of these sources, which had enough data to perform these estimates. The derived $e$-folding timescales were remarkably short, only of the order of days. The average radio luminosities based on the 37 GHz data are of order $10^{42}$–$10^{43}$ erg s$^{-1}$.

Järvelä et al. [35] discuss several possible origins for this extreme behaviour, and after excluding most of them, three alternatives could not be discarded: i) jet or jet base interacting with a star or a BLR cloud; ii) relativistic jets in the presence of extreme free-free absorption; and iii) magnetic reconnection. All of these exhibit distinct multiwavelength characteristics, such as timescales and correlations, particularly in radio and X-rays. Thus, simultaneous radio and X-ray observations could help shed light on the origin of this extreme phenomenon.

To that end, in this study, we analyse simultaneous radio and X-ray observations, complemented by ultraviolet and optical data. In late May 2023, J1228+5017 exhibited a high flux density flare at 37 GHz (780$\pm$120 mJy). Due to this, we triggered the JVLA observations and successfully proposed Target-of-Opportunity observations with the Neil Gehrels Swift Observatory ([20]) to follow up on this flare. OVRO was already monitoring this source. With these data, our goal is to investigate the multiwavelength behaviour of this source; for example, we are interested in whether these bands show any correlated variability, similar to what was previously observed in one of these sources ([70]), or if the higher-energy optical/ultraviolet/X-ray bands show any signs of significant variability. Furthermore, we aim to further constrain the multiband timescales associated with the flaring, since those might help us narrow down the possible explanations for this behaviour. In Sect. 2 we introduce the source, in Sect. 3 we describe the observations and the data reduction and analysis, in Sect. 4 we present our results, which are discussed in Sect. 5. Finally, we summarise our results and present our conclusions in Sect. 6. Throughout this paper, we adopt a standard $\Lambda$CDM cosmology, with a Hubble constant $H_{0} = 72$ km s$^{-1}$ Mpc$^{-1}$, and $\Omega_{\Lambda} = 0.73$ ([75]). For spectral indices, $\alpha$, we adopt the $S_{\nu} \propto\nu^{\alpha}$ convention, where $S_{\nu}$ is the flux density at frequency $\nu$.

J1228+5017

J1228+5017 (RA 12:28:44.82, Dec +50:17:51.24) is an NLS1 galaxy at a redshift of $z=0.2627$². It is hosted by a disk-like galaxy ([82]) and was chosen for the MRO monitoring due to its very dense large-scale (Mpc-scale) environment ([34, 44]). It was first detected at 37 GHz in 2015, and before the 2023 flare, the last detection was in 2019 [35]. OVRO has also detected it at 15 GHz both before and after the 2023 flare (Crepaldi et al. in prep.). It was also identified as a possible gamma-ray emitter by Lähteenmäki et al. [44], who reported a test statistic value of 20, corresponding to a 4.5$\sigma$ tentative detection, yet to be confirmed.

Crepaldi et al. [16] obtained and modelled a very high-quality (signal-to-noise $\sim65$ in the continuum) optical spectrum of J1228+5017. They estimated its black hole mass to be $(9.21 \pm4.96) \times10^{6} M_{\odot}$, and the Eddington ratio to be $0.11 \pm0.02$ or $0.94 \pm0.19$, depending on which method was used to estimate the continuum luminosity. The smaller value results from using the H$\beta_{\mathrm{broad}}$ luminosity to estimate the monochromatic luminosity at 5100 Å, and the larger value from directly measuring the 5100 Å luminosity. Given its extremely strong Fe II emission, with an R4570 parameter (strength of Fe II relative to total H$\beta$) of 2.81, the latter Eddington ratio estimate is more likely to be closer to the true value, as R4570 is known to correlate with the Eddington ratio [49]. This interpretation is supported by the strong outflows seen in the [O III]$\lambda\lambda4959,5007$ lines, often driven by the strong radiation pressure. Modelling these lines requires two wing components, which are blueshifted by $-163$ and $-371$ km s$^{-1}$, indicating strong outflows (see Figure 2 of Crepaldi et al. [16]).

There is only one archival X-ray detection of J1228+5017 from ROSAT at 0.0452 counts per second (cps). Thus, we had no prior knowledge of its X-ray properties. However, all in all, J1228+5017 is a stereotypical NLS1 with a low-mass supermassive black hole, high Eddington ratio, strong Fe II emission, and [O III] emission lines with strong outflows.

Data

We have simultaneous data from two single-dish telescopes, the JVLA, and the X-ray Telescope (XRT, Burrows et al. [12]) and Ultra-Violet/Optical Telescope (UVOT, Roming et al. [71]) onboard Swift. The single-dish radio telescopes MRO and OVRO operate at 37 and 15 GHz, respectively. Our JVLA data are in X, K, and Ka bands, centred at 10, 22, and 33 GHz, respectively. XRT data cover the 0.3-10 keV energy band, and our UVOT observations are in V, B, U, W1, M2, and W2 bands. The data included in this work were collected between May 24 and July 13, 2023.

Metsähovi Radio Observatory

The measurements included in this study are part of the large ongoing AGN monitoring programme at 22 and 37 GHz with the 13.7-m radio telescope at MRO. The observations are made with a 1 GHz-band dual-beam receiver centred at 36.8 GHz. The beam full-width at half power is 144 arcsec. The observations are on–on observations, alternating the source and the sky in each feed horn. A typical integration time to obtain one flux density data point of a faint source is 1800 s. The sensitivity is limited by sky noise due to the location of the telescope, and it has been experimentally shown that the results do not significantly improve after the maximum integration time of 1800 s. The detection limit of the telescope at 37 GHz is of the order of 200 mJy under optimal conditions. Data points with an S/N $< 4$ are handled as non-detections. The flux density scale is set by observations of DR 21. Sources NGC 7027, 3C 274, and 3C 84 are used as secondary calibrators. A detailed description of the data reduction and analysis is given in Teraesranta et al. [78]. The error estimate for the flux density includes contributions from the measurement rms and the uncertainty in the absolute calibration. MRO detections for J1228+5017 are shown in Table 1 and the MRO light curve, including upper limits, is shown in Figure 2.

date (dec. year)	date and time (yyyy-mm-dd hh:mm:ss)	$S_{\mathrm{int}}$ (mJy)
2023.405411	2023-05-28 23:23:53	$780 \pm120$
2023.427261	2023-06-05 22:48:19	$440 \pm80$
2023.432697	2023-06-07 22:25:21	$500 \pm110$
2023.442703	2023-06-11 14:04:43	$870 \pm170$
2023.443627	2023-06-11 22:10:18	$610 \pm120$
2023.450887	2023-06-14 13:46:12	$630 \pm150$
2023.454555	2023-06-15 21:53:54	$600 \pm130$
2023.476434	2023-06-23 21:33:45	$490 \pm90$
2023.506495	2023-07-04 20:54:00	$480 \pm90$

Owens Valley Radio Observatory

The 15 GHz observations were carried out as part of the general OVRO 40 m radio telescope AGN monitoring programme. This telescope uses off-axis dual-beam optics and a cryogenic receiver with a 15.0 GHz centre frequency and 3 GHz bandwidth. The beam full-width at half power is 157 arcsec. The observations are carried out in an on–on fashion to remove atmospheric and ground contamination. In May 2014, a new pseudo-correlation receiver was installed on the 40 m telescope, and the fast gain variations are corrected using a 180-degree phase switch. Calibration is achieved using a temperature-stable diode noise source to remove receiver gain drifts, and the flux density scale is derived from observations of 3C 286, assuming the Baars et al. [3] value of 3.44 Jy at 15.0 GHz. The systematic uncertainty of about 5 per cent in the flux density scale is included in the error bars. The upper limits are $4\sigma$ upper limits based on the measurement rms. Complete details of the reduction and calibration procedure are found in Richards et al. [69], and more details specific to the NLS1 observations are given in Appendix A3 Järvelä et al. [35]. J1228+5017 has been detected by OVRO, but we did not detect it during the period used in this paper. The OVRO upper limits are shown in Figure 2.

Karl G. Jansky Very Large Array

We observed J1228+5017 with the JVLA a total of seven times in three different bands, X, K, and Ka, centred at 10, 22, and 33 GHz, respectively (Project VLA/23A-043, PI Järvelä). The dates and times of the JVLA observations are given in Table 2. The array configuration changed over our observing period: for the first three epochs, the array was in BnA, the next three in BnA$\rightarrow$A, and the last observation in A configuration. The resolutions of these configurations are comparable, and the source was point-like in most observations; thus, we have no reason to assume that configuration changes would have significantly affected the measured flux densities. The total bandwidth was 4 GHz in X, and 8 GHz in K and Ka bands, each band divided into 128 MHz subbands, consisting of 64 channels of 2 MHz. The NLS1 3C 286 (Berton et al. [7]) was used as the bandpass and flux density calibrator for each source, and each source had an individual nearby bright source used as the complex gain calibrator. The pointing offset calibration was done either at 3C 286 or the current complex gain calibrator. The expected thermal noise levels were 13, 18, and 22 $\mu$Jy in X, K, and Ka bands, respectively. In most cases, we reached even lower rms levels (see Table 2).

Table 2: Interferometric data for J1228 in X, K and Ka bands. Columns: (1) date in decimal years; (2) date and the observation end time; (3) peak flux density and its error (from the 2D Gaussian fit), not including the flux density calibration uncertainty; (4) integrated flux density and its error; (5) rms level of the observation; (6) clean beam size; (7) beam position angle; (8) array configuration.

X
date (dec. year)	date and time (yyyy-mm-dd hh:mm:ss)	$S_{\mathrm{peak}}$ ($\mu$Jy beam$^{-1}$)	$S_{\mathrm{int}}$ ($\mu$Jy)	rms ($\mu$Jy beam$^{-1}$)	beam size (arcsec)	beam PA (deg)	array configuration
2023.449861 (VLA 1)	2023-06-14 04:47:11	163 ± 10		12	0.64 × 0.49	67.8	BnA
2023.449970 (VLA 2)	2023-06-14 05:44:08	176 ± 10		12	0.74 × 0.43	60.0	BnA
2023.452638 (VLA 3)	2023-06-15 05:06:30	158 ± 9		11	0.66 × 0.48	66.0	BnA
2023.471926 (VLA 4)	2023-06-22 06:04:25	182 ± 11		13	0.94 × 0.26	52.2	BnA→A
2023.488120 (VLA 5)	2023-06-28 03:56:06	163 ± 9		11	0.26 × 0.22	−64.3	BnA→A
2023.488229 (VLA 6)	2023-06-28 04:53:02	152 ± 7		13	0.28 × 0.21	−76.3	BnA→A
2023.496394 (VLA 7)	2023-07-01 04:24:50	153 ± 9		11	0.26 × 0.22	−71.8	A
K
date (dec. year)	date and time (yyyy-mm-dd hh:mm:ss)	$S_{\mathrm{peak}}$ ($\mu$Jy beam$^{-1}$)	$S_{\mathrm{int}}$ ($\mu$Jy)	rms ($\mu$Jy beam$^{-1}$)	beam size (arcsec)	beam PA (deg)	array configuration
2023.449861 (VLA 1)	2023-06-14 04:47:11	120 ± 7		10	0.29 × 0.22	71.5	BnA
2023.449970 (VLA 2)	2023-06-14 05:44:08	168 ± 8		11	0.30 × 0.21	71.1	BnA
2023.452638 (VLA 3)	2023-06-15 05:06:30	119 ± 7		10	0.30 × 0.24	63.1	BnA
2023.471926 (VLA 4)	2023-06-22 06:04:25	128 ± 12		15	0.38 × 0.12	54.0	BnA→A
2023.488120 (VLA 5)	2023-06-28 03:56:06	116 ± 9		11	0.12 × 0.09	−64.8	BnA→A
2023.488229 (VLA 6)	2023-06-28 04:53:02	80 ± 11		13	0.13 × 0.10	−77.7	BnA→A
2023.496394 (VLA 7)	2023-07-01 04:24:50	130 ± 10		10	0.12 × 0.09	−68.4	A
Ka
date (dec. year)	date and time (yyyy-mm-dd hh:mm:ss)	$S_{\mathrm{peak}}$ ($\mu$Jy beam$^{-1}$)	$S_{\mathrm{int}}$ ($\mu$Jy)	rms ($\mu$Jy beam$^{-1}$)	beam size (arcsec)	beam PA (deg)	array configuration
2023.449861 (VLA 1)	2023-06-14 04:47:11	161 ± 11		14	0.19 × 0.16	79.1	BnA
2023.449970 (VLA 2)	2023-06-14 05:44:08	190 ± 13		14	0.20 × 0.14	70.6	BnA
2023.452638 (VLA 3)	2023-06-15 05:06:30	110 ± 11		13	0.20 × 0.17	61.7	BnA
2023.471926 (VLA 4)	2023-06-22 06:04:25	135 ± 13		16	0.26 × 0.08	53.8	BnA→A
2023.488120 (VLA 5)	2023-06-28 03:56:06	96 ± 11		15	0.08 × 0.06	−61.7	BnA→A
2023.488229 (VLA 6)	2023-06-28 04:53:02	121 ± 11	130 ± 20	15	0.08 × 0.06	−77.9	BnA→A
2023.496394 (VLA 7)	2023-07-01 04:24:50	100 ± 11	120 ± 18	13	0.09 × 0.06	−67.2	A

We used the Science Ready Data Products (SRDP) provided by the NRAO. The data were calibrated using the VLA Imaging Pipeline 2022.2.0.64. After this, the data were checked manually, and any remaining bad data were flagged, producing the SRDP measurement set for each source. We also performed a manual check of the data, but no additional flagging was required. In the subsequent steps, we used the Common Astronomy Software Applications (CASA) version 6.4.1.12. We averaged the J1228+5017 data in time (timebin = 10 s) and in frequency (width = 64, averaging 64 channels to form one output channel per spectral window).

The radio maps of our sources were produced using the tclean algorithm with interactive cleaning. The cell size was chosen to ensure the synthesised beam is properly sampled; we used cell sizes of 0.025, 0.010, and 0.007 arcsec in the X, K, and Ka bands, respectively. The image size in the X band was 512×512 pixels, whereas we used 1024×1024 in the K and Ka bands. We used Briggs weighting, with robust = 1.8, in all cases. We used the mtmfs deconvolver with nterms = 2 and scales = 0 in case J1228+5017 would be bright enough to produce spatially resolved in-band spectral index maps, which turned out not to be the case. We applied the wide-band primary beam correction separately using widebandpbcor.

We fitted the source in each map using a 2D Gaussian to obtain the central coordinates, the peak flux density, and its error. In the case of extended emission (only two Ka band observations), we measured the emission within the 3σ contour and estimated its error by multiplying the rms by the square root of the emitting region area, expressed in beams. The rms for each map was measured in an empty region of sky far from the central source. Details of the measurements are given in Table 2, the radio spectra for the seven epochs are shown in Figure 3, and the JVLA light curves are shown in Figure 2. Since the JVLA flux density scale calibration accuracy is of order 5% in X band, and 10-15% in K and Ka bands, we took this into account in the error bars shown in Figure 3 by adding in quadrature. The errors reported in Table 2 as well as the errors shown in Figure 2 do not include this additional source of error.

Neil Gehrels Swift Observatory

The Swift data were collected through target of opportunity observations (Target ID 16062, PI: P. Romano, 13 observations) from 2023-05-30 to 2023-07-04 with XRT (Burrows et al. [12]) and the UVOT (Roming et al. [71], with the “U+B+V+all U” filter setup), for a total exposure of 34.7 ks.

The Swift/XRT data were reduced and analysed by using standard procedures within FTOOLS³ (v6.36), and the matching calibration database CALDB⁴. In particular, the XRT data were processed and filtered with the task xrtpipeline (v0.13.7). The XRT light curve was generated in the 0.3–10 keV energy range with one point per observation. We used sosta within XIMAGE with a background measured in an annular region with an inner radius of 60 pixels and an outer radius of 110 pixels centred at the source position. The appropriate exposure maps were used to correct for PSF losses and vignetting. We also extracted a total spectrum by adding all XRT observations (34.7 ks). For the spectroscopic analysis, source events were extracted from a circular region with a 15-pixel radius; ancillary response files were generated with the task xrtmkarf to account for different extraction regions, vignetting, and PSF corrections; the spectral redistribution matrices in CALDB were used. We fit the total spectrum with XSPEC (v. 12.15.1, Arnaud [2]) with a TBabszTBabsPOWERLAW, that is a simple power law model corrected for an equivalent hydrogen column corresponding to the Galactic value (fixed $N_{\rm H}^{\rm Gal} = 1.56 \times10^{20}\ {\rm cm}^{-2}$, HI4PI Collaboration et al. [27]) and with the TBABS (Wilms et al. [87]) model and a further absorption column at fixed redshift $z = 0.262$. We obtain a photon index $\Gamma= 2.33^{+0.20}_{-0.18}$ and an observed (unabsorbed) 0.3-10 keV flux of $(6.6 \pm0.6) \times10^{-13}$ $(7.5^{+0.9}_{-0.7} \times10^{-13})$ erg cm$^{-2}$ s$^{-1}$ (statistics/dof = 156.36/234).

The analysis of the UVOT data, collected in all filters (optical $v$, $b$, $u$ and UV $w1$, $m2$, $w2$), was performed with the tasks UVOTMAGHIST, UVOTIMSUM, and UVOTSOURCE within FTOOLS (v6.34), which offer a check for frames affected by small scale sensitivity, SSS, issues⁵. We used UVOTMAGHIST to identify the frames of each UVOT image affected by SSS issues and to generate a light curve in each filter at the frame level (the best time resolution available), UVOTIMSUM to sum the sky images for UVOTSOURCE to calculate the magnitude of the source through aperture photometry within a circular region centred on the best source position by applying the required corrections related to the specific detector characteristics. We adopted circular regions with radii of 5″ and 20″ for the source and background, respectively.

The Swift/XRT observing epochs and exposure times are shown in Table A.1, and the light curves for XRT and UVOT can be found in Figure 2.

Results

Note on measured variability

Before presenting the results, we briefly discuss the variability measures and statistics used. We performed the weighted $\chi^2_{\rm r}$ test, as well as calculated the $C$ and $F$ parameters for all light curves (Singh & Meintjes [73]). These tests indicate that only the XRT and W2 light curves are variable. However, the number of data points in the light curves ranges from seven to 13, leading to low degrees of freedom. This directly determines the $\chi^2$ distribution and adversely impacts the power and validity of these tests. Moreover, 7-13 observations cannot properly represent the underlying flux density distribution, and generally it is not advised to use statistical tests based on the normal distribution of the data with fewer than 30 data points. With that being said, based on these tests, we cannot rule out the possibility that all the light curves, except for XRT and W2, are constant. However, several light curves do appear variable, as we will discuss in the following sections.

Single-dish observations

All the follow-up observations were triggered by an MRO detection on May 28, 2023. The flux density of this detection was 780 ± 120 mJy. During the observing period shown in Figure 2, J1228+5017 was detected a total of nine times, with $S_{\nu,\max} = 870 \pm170$ mJy, $S_{\nu,\min} = 440 \pm80$ mJy and $S_{\nu,\mathrm{ave}} = 600$ mJy. During this period, J1228+5017 was observed 61 times, yielding a detection rate of 15%. Upper limits were derived for the 52 non-detections. The MRO light curve shows repeating flaring behaviour at 37 GHz; on average, the upper limits are lower than the detections, suggesting the source would have been detected had it maintained its high flux density. However, based on the 37 GHz observations, we cannot say much about the flux density levels between the detections; in principle, the source could have been just below the detection limit.

Interestingly, OVRO did not detect the source during this period. The 15 GHz observations yielded only five upper limits, on the order of a few tens of mJy. On June 15, an OVRO non-detection was followed by an MRO detection only 1.2 hours later, indicating a very short timescale, assuming that this phenomenon is broadband. The OVRO beam size is close to that of

MRO, so OVRO upper limits can be used to estimate the low flux density baseline for the MRO flux densities by assuming some reasonable spectral index between the two frequencies. These will be discussed more in Sect. 4.5.

Interferometric observations

Our arcsec-resolution radio data consist of seven epochs of JVLA observations in the X, K, and Ka bands from June 14 to July 1, 2023. J1228+5017 was detected at $> 6\sigma$ level in all bands at all epochs (see Table 2). It was otherwise consistent with a point-like morphology except for the last two Ka band observations, where it was marginally resolved. However, within the errors, its integrated flux density remains consistent with the corresponding peak flux density, indicating that the amount of any resolved emission is minimal. JVLA also changed its configuration from BnA to A during these observations, with the last observations having the highest resolution. There are no trends of decreasing flux density moving from BnA to A configuration, ruling out resolved-out emission having a significant impact on the observed variability.

All bands seem to show variability over the observing period, even if, based on the variability test, a constant flux density cannot be ruled out. The X band shows the least change in flux density, but even in this case, the minimum and maximum flux densities differ, even when accounting for the flux density scale accuracy, as can be seen in Figure 1. The variability is more pronounced in K and Ka bands, with the maximum flux density being approximately double the minimum flux density. As mentioned earlier, neither the array configuration nor the clean beam size appears to correlate significantly with flux density.

The average peak flux densities in X, K, and Ka bands are $164 \pm4$, $123 \pm10$, and $130 \pm13$ $\mu$Jy beam$^{-1}$, respectively. We estimate the error as the standard error of the mean. Compared to the peak flux densities reported by Järvelä et al. [35], the X and Ka band values are statistically higher; they report X, K, and Ka band peak flux densities of $128 \pm5$, $120 \pm7$, and $102 \pm11$ $\mu$Jy beam$^{-1}$.

As a result of possible variability in the K and Ka bands, the spectral shape in the 10-33 GHz range also appears to change, showing flat, concave, and steep shapes. Still, even at its steepest at epoch 5, the spectral index is $\alpha=-0.44$, and thus the spectrum is formally flat.⁶ This is consistent with the earlier broadband spectrum shown in Figure 1, which indicates that the spectrum flattens above 10 GHz.

There are two instances where we have detection at 37 GHz from MRO and in Ka band from JVLA within 12 hours of each other. These consecutive observations will be discussed in Sect. 4.5.

Optical, ultraviolet, and X-ray observations

We have 13 epochs of Swift data obtained with XRT and UVOT. The XRT energy range is 0.3–10 keV, and the UVOT bands used are V, B, U, W1, M2, and W2, with central wavelengths of 5468, 4392, 3465, 2600, 2246, and 1928 Å, respectively. The light curves are shown in Figure 2.

The optical and UV light curves seem to show variability by a factor of 1.5–2, but we identify no obvious flare-like events. V and B bands seem to be the most variable, which might be due to the fact that the Balmer discontinuity falls in the B band, and the H$\gamma$ and approximately half of the blueward Fe II bump fall in the V band. However, based on the archival optical spectra of J1228+5017, its spectrum has stayed remarkably similar over the past 20 years. The only statistically variable optical/UV light curve is the W2 band; this is likely due to its very small errors, as variability tests are highly sensitive to the size of the measurement errors.

The X-ray light curve is slightly more variable, also according to the $\chi^2$ variance test, with variability by a factor of $\sim$3. Similar to the V and B band light curves, the X-rays seem to be more elevated during the first half of the observing period, but no flare-like events with unusually high amplitude can be seen. We observe no clear, similar trends when comparing the optical, UV, or X-ray light curves to the radio light curves.

Variability timescales

We can use consecutive detections, or a nearby detection and an upper limit, to estimate variability properties, such as timescales, variability brightness temperatures, and variability Doppler factors of the flares (Valtaoja et al. 1999; Hovatta et al. [30]). In the following, we assume that the flux density differences between bands and/or instruments trace the evolution of the same variable component. We will consider the following cases: i) MRO 37 GHz detection and JVLA Ka band detection within 24 hours of each other; and ii) MRO detection and OVRO non-detection within 24 hours of each other. We will also consider two consecutive JVLA or MRO detections, but in those cases, the flux densities are similar when accounting for the errors. Before performing any calculations, we make some remarks on these cases and explain the assumptions we will use.

When using MRO and JVLA data, the resolutions differ significantly: 144 arcsec for MRO and much less than 1 arcsec for the JVLA in the Ka band (33 GHz). The largest angular scale detectable with the JVLA in A configuration in Ka band is 1.6 arcsec, which translates to 6.3 kpc at the redshift of our source. We do not know how much of the emission is resolved out with the JVLA, which is why, in these cases, we will use the closest OVRO upper limit and a spectral index of $\alpha= 0$ to estimate the baseline flux density at 37 GHz. This approach is based on the similar beam sizes of MRO and OVRO — they should probe approximately the same radio-emitting regions, with the flux density modulated by the spectral index — and on the fact that, due to its spatial extent, the resolved-out emission cannot vary on short timescales. We adopt a spectral index of 0, which is close to the spectral slope of free-free emission (i.e. −0.1). As far as we know, there are no emission mechanisms that produce inverted spectra at kpc-scales. For JVLA flux density errors, we will use values that account for the flux density scale uncertainty.

Since MRO and OVRO are observing at different frequencies, when comparing an MRO detection and an OVRO upper limit, we will extrapolate the OVRO upper limit to 37 GHz using an extreme spectral index of $\alpha= 4$. We adopt this value to represent the worst-case scenario, for example, extreme synchrotron self-absorption coupled with external free-free absorption. In this case, we assume the flares are broadband, which may not be true. We also want to remind the reader that upper limits are not detections and should not be treated as such; the results using an upper limit and extrapolation of data will be less accurate and reliable than the ones derived with detections.

With these caveats and assumptions in mind, we can estimate the timescales, variability brightness temperatures, and variability Doppler factors of the flares. For simplicity, we used the same equation for both rising and decaying flares:

where $\Delta S_{\max}$ is the maximum amplitude of the flare in Jy, after subtracting the baseline flux density level, $S_{\mathrm{b}}$, $t_{\max}$ is the epoch of the peak of the flare, and $\tau$ is the rise or decay time of the flare expressed in days (e-folding timescale). In our case, we only have two measurements per flare, so we can simplify $|t-t_{\max}|$ to $\Delta t$, which is the time between the middle points of the observations.

The variability brightness temperature, $T_{\mathrm{b,var}}$, in the source proper frame was estimated with:

where $\nu$ is the observed frequency in GHz, $d_{\mathrm{L}}$ is the luminosity distance in metres, and $\Delta S_{\max}$ and $\tau$ are defined in Eq. (1). The numerical factor corresponds to using $H_0 = 72\ \mathrm{km\ s^{-1}\ Mpc^{-1}}$, and $\Omega_\Lambda= 0.73$, and assuming that the source is a homogeneous sphere. Estimating the brightness temperature from flux density variability relies on a causality argument, so these values are lower limits. We calculated the variability brightness temperatures for all flares with $\tau$ values. We note that the brightness temperatures derived from variability are systematically larger by a factor of $\delta^3$, where $\delta$ is the Doppler factor, than those obtained directly from high-resolution measurements due to the different dependence on the Doppler factor.

Finally, we can use the variability brightness temperature to estimate the variability Doppler factor, assuming we know the intrinsic brightness temperature, $T_{\mathrm{b,int}}$:

For the intrinsic brightness temperature, we use $5 \times10^{10}\ \mathrm{K}$ [68, 45] (Homan et al. 2021), which assumes equipartition between the energy densities of the magnetic field and the radiating particles. We do not know whether J1228+5017 is in equipartition and therefore cannot assess the accuracy of the Doppler factor estimates.

The $e$-folding timescales, variability brightness temperatures and variability Doppler factors for three flares are shown in Table 3. In two cases, J1228+5017 was detected in Ka band (33 GHz) and shortly after at 37 GHz. In one case, a 37 GHz detection was followed by a 15 GHz non-detection only 1.2 h later. The $e$-folding timescales in all cases are on the order of hours, even when accounting for errors. The derived variability brightness temperatures are of the order of $10^{17} - 10^{18}\ \mathrm{K}$. These values far exceed the Compton catastrophe threshold of $10^{12}\ \mathrm{K}$, but can be partly explained by extreme Doppler factors, due to the dependence of the $\delta^3$ term. Resulting from the high variability brightness temperatures, the derived variability Doppler factors are of the order of a few hundred.

facilities	date	$\Delta t$	$S_{\mathrm{begin}}$	$S_{\mathrm{end}}$	$\tau$	$T_{\mathrm{b,var}}$	$\delta_{\mathrm{var}}$
	(yyyy-mm-dd)	(hours)	(mJy)	(mJy)	(hours)	($10^{18}$ K)
VLA Ka – MRO	2023-06-14	8.9672	$0.19 \pm0.03^{a}$	$618.00^{b} \pm150.00$	$1.1088^{+0.0615}_{-0.0509}$	$4.20^{+1.53}_{-1.35}$	$438^{+48}_{-53}$
VLA Ka – MRO	2023-06-15	17.72278	$0.11 \pm0.02^{a}$	$588.00^{b} \pm130.00$	$2.0646^{+0.1054}_{-0.0920}$	$1.15^{+0.40}_{-0.34}$	$284^{+30}_{-31}$
MRO – OVRO	2023-06-15	1.1919	$588.00^{b} \pm130.00$	$424.00^{c}$	$3.6450^{+11.8070}_{-1.3822}$	$0.37^{+0.80}_{-0.35}$	$195^{+91}_{-121}$

Using a simple light-travel-time argument to estimate the size of the radio-emitting region, its size needs to be $r < c\tau\delta/(1+z)$. For the shortest $e$-folding timescale, this gives $r < 0.03\ \mathrm{mpc} \times\delta$. If the Doppler factor exceeds one hundred, the emitting region would be on the order of a few mpc, or 3-4 light-days; comparable to the size of the accretion disk. Adopting the Crepaldi et al. [16] value for the black hole mass, $9.21 \times10^6 M_\odot$, the gravitational radius, $r_{\mathrm{g}}$, of J1228+5017 is $4.4 \times10^{-4}\ \mathrm{mpc}$. It should be remembered that the derived timescales are upper limits, and the actual variability might occur faster, indicating smaller emitting regions.

Discussion

A thorough discussion of alternative explanations, intrinsic and external to the source, is given by Järvelä et al. [35] and will not be repeated here. Instead, we focus on interpreting the multi-wavelength monitoring of J1228+5017 and the extremely short timescales, coupled with high-amplitude variability, $T_{\mathrm{B,var}}$ and $\delta_{\mathrm{var}}$.

Silence of the JVLA

The JVLA observations pose two interesting questions. First, why are these flares not detected by the JVLA, and second, what is the origin of the high-frequency flat-spectrum radio emission seen by the JVLA?

The JVLA observed J1228+5017 seven times over approximately two weeks. It coincides with a period when MRO observed the source 28 times and detected it three times, giving a detection fraction of $\sim10\%$ (see Järvelä et al. [35] for a thorough discussion about the reliability of the MRO detections). With a $10\%$ duty cycle, it is not an unreasonable assumption that the seven JVLA observations just missed the flares. However, given MRO’s moderate sensitivity, we would expect a higher duty cycle when observed with the higher-sensitivity JVLA, since it should also detect rising or decaying flares that remain below MRO’s detection threshold. Since this is not the case, it may indicate even shorter timescales or flare profiles not well represented by exponential rise and decay. The nature of the flares is discussed in more detail in Sect. 5.3.

Given the low flux densities between 10 and 33 GHz, it is worth exploring whether star formation, rather than an AGN, could be responsible for the radio emission. At these frequencies, arcsec-scale AGN radio spectra are often dominated by a steep ($\alpha\approx-0.7$) optically thin synchrotron component (Panessa et al. [63]), whereas in the presence of star formation, toward high frequencies, the flat free-free emission ($\alpha\approx-0.1$) starts to flatten the steep supernova remnant synchrotron spectrum ($\alpha\approx-0.8$ Klein et al. [39]). This spectrum often also has a break around 1-12 GHz, further decreasing its contribution at higher frequencies (Klein et al. [39]). At 33 GHz, the thermal fraction is already around 80% (Digman et al. 2025), indicating that the radio spectrum should be flat. This alone does not prove that star formation is present; thus, we will utilise additional star formation proxies to estimate its contribution. First, we can estimate the star formation rate (SFR) required to produce the 33 GHz radio emission following Murphy et al. [53]. We assume $T_{\mathrm{e}} = 10^4\ \mathrm{K}$ and a non-thermal spectral index, $\alpha_{\mathrm{NH}} = -0.7$. We calculate the luminosity, $L_\nu$, from the average flux density of the last three epochs ($105.7 \pm19.1\ \mathrm{\mu Jy}$), since these observations have the smallest clean beam, thus providing the best constraints on the size of the emitting region. Based on the resulting luminosity, $L_{33\ \mathrm{GHz}} = (1.70 \pm0.3) \times10^{29}\ \mathrm{erg\ s^{-1}\ Hz^{-1}}$, the SFR would be $\sim60 \pm10\ M_\odot\ \mathrm{yr^{-1}}$ within the innermost $(0.36 \times0.24)\ \mathrm{kpc}^2$. This would yield an SFR surface density of $\sim694\ M_\odot\ \mathrm{yr^{-1}\ kpc^{-2}}$, which is unreasonably high compared to the average values in the local Universe, which are around $0.5\ M_\odot\ \mathrm{yr^{-1}\ kpc^{-2}}$ (Winkel et al. [88]), and thus it is unlikely that star formation has a significant contribution to the radio emission at 33 GHz.

We can also use Eq. 13 by Murphy et al. [53] to estimate the SFR from mid-infrared luminosity; instead of 24 $\mathrm{\mu m}$, we will use the 22 $\mathrm{\mu m}$ value from the Wide-field Infrared Survey Explorer (WISE). However, it is worth noting that the AGN also contributes to the mid-infrared luminosity, and when polar dust is present, it can produce a mid-infrared excess, mimicking observational signatures of star formation (Järvelä et al. [33]). Keeping that in mind, interestingly, estimating the SFR from the 22 $\mathrm{\mu m}$ luminosity gives exactly the same value of $\sim60\ M_\odot\ \mathrm{yr^{-1}}$, but in this case, due to the resolution of WISE, for the whole galaxy. This is not an unreasonable value, especially for NLS1s, which are known to exhibit enhanced star formation (Sani et al. [72]). Alternative ways to assess the strength of star formation are the infrared colours and the $q22$ parameter. A very red W3-W4 colour (12$\mathrm{\mu m}$–22$\mathrm{\mu m}$), exceeding 2.5, indicates a strong star formation contribution to the infrared emission (Caccianiga et al. [13]), since the emission originating from star formation increases faster relative to the emission of AGN origin. Another parameter they suggest is $q22$, defined as $\log(S_{22\ \mathrm{\mu m}} / S_{1.4\ \mathrm{GHz}})$. This parameter is analogous to the classical far-infrared - radio correlation parameter, but uses the WISE $22\ \mu\mathrm{m}$ band instead of far-infrared measurements. According to them, $q_{22}$ exceeding unity can be a sign of enhanced star formation. For J1228+5017, W3-W4 = 2.45, and $q_{22} = 1.9$, indicating possibly enhanced star formation in the host.

All these proxies use the $22\ \mu\mathrm{m}$ measurement, so the question of what fraction of this emission comes from the AGN versus star formation remains. As mentioned earlier, AGN polar dust can enhance mid-infrared emission, and especially so in sources with strong narrow-line region outflows [48]. Indeed, J1228+5017’s [O III]$\lambda\lambda$4959,5007 emission lines exhibit a strong, high-velocity outflow (see Figure 2 in Crepaldi et al. [16]). Two wing components are required to properly model the outflow, with the higher-velocity component showing a $(-371 \pm35)\ \mathrm{km\ s^{-1}}$ blueshift. The main mechanisms to drive these outflows are the accretion disk winds [66] and (relativistic) jets [77, 55, 8]. In the case of J1228+5017, no persistent jets seem to be present, so the NLR outflow is likely driven by radiation pressure induced by the high Eddington ratio. All things considered, the current data, unfortunately, do not allow us to determine the origin of the $22\ \mu\mathrm{m}$ emission.

However, other properties of our JVLA data do not indicate a star formation origin for the high-frequency radio emission. As can be seen in the light curves in Figure 2 and the spectra in Figure 1, the K and Ka band emission seems to be variable. As we do not expect to see variability from star formation at these timescales, it is safe to assume that at least the varying fraction of the flux density is of AGN origin. Thus, if star formation is not likely to dominate the highest observed frequencies, we need to turn our attention back to AGN-produced emission. In the absence of extended optically thin emission, detecting the optically thick jet base can yield a flat spectrum [63]. However, we would have expected to detect it also in the mas-resolution VLBA observations at 15 GHz. The 1$\sigma$ rms of those observations was 58 $\mu$Jy, and thus any possible jet base emission was probably under the detection threshold [35].

Another component that should start to become detectable toward higher radio frequencies is the tentative radio corona [46]. Its possible physical connection to the jet base is still unclear, but for clarity, we will discuss it as a separate component. The radio corona is thought to be magnetically heated and, at its smallest scales at hundreds of GHz, cospatial with the X-ray corona, with a size around some tens of $r_{\mathrm{g}}$ — similarly to coronally active stars. These stars, as well as several non-jetted AGN, follow the Güdel-Benz relation: $L_{\mathrm{R}}/L_{\mathrm{X}} \sim10^{-5}$ [25]. We estimate the X-ray luminosity from the Swift observations assuming 0.03 cps between 0.3–10 keV and a photon index of 2.5. For radio, we use the earlier calculated $L_{33\ \mathrm{GHz}}$ to calculate $\nu L_{\nu}$. This yields $L_{\mathrm{R}}/L_{\mathrm{X}} = 4 \times10^{-5}$, being consistent with the Güdel-Benz relation. Assuming equipartition and a power-law electron energy distribution with $p=2$, we can estimate the size of the radio-emitting region using Eq. 22 by Laor & Behar [46]:

where $R_{\mathrm{pc}}^{\mathrm{RS}}$ is the size of the radio sphere in pc, $L_{30}$ is $L_{\nu}/10^{30}\ \mathrm{erg\ s^{-1}\ Hz^{-1}}$, $L_{46}$ is $L_{\mathrm{bol}}/10^{46}$ and $\nu_{\mathrm{GHz}}$ is the observed frequency in GHz. To estimate $L_{30}$, we use $L_{33\ \mathrm{GHz}}$, and for the bolometric luminosity we use the value derived from the monochromatic 5100 Å luminosity by Crepaldi et al. [16]. J1228+5017 does not appear to host a persistent jet, and based on the multiwavelength observations so far, the extreme radio variability does not seem to have an optical counterpart; therefore, we assume that the 5100 Å luminosity is not contaminated by this activity. With these values, we obtain $R_{\mathrm{pc}}^{\mathrm{RS}} = 0.006$ pc or 6 mpc, which corresponds to $\sim1.4 \times10^{4}\ r_{\mathrm{g}}$, or around 7 light-days. Taken at face value, this size is likely to be too large to be related to the X-ray or the radio corona. However, it is consistent with the size of the emitting region derived from the variability timescale, $r < 0.03 \times\delta$, assuming $\delta> 100$, but it should be noted that these estimates are based on very different assumptions. Both assume equipartition, but flare timescales assume relativistically boosted emission, whereas the coronal emission is expected to be close to isotropic. This could be a coincidence, or it could point to the low- and flaring-state emission originating from the same region.

Lack of multiwavelength trends

Independently of the origin of the radio emission, the jet (base) or the corona, in several cases, we would expect to see correlated multiwavelength variability induced by the couplings between the accretion disk, the X-ray corona, and the radio corona or the jet. In another similar source, SDSS J164100.10+345452.7, X-rays brightened after a 37 GHz flare, possibly suggesting a connection between the two bands [70]. However, in the case of J1228+5017, based on the multiwavelength light curves in Figure 2, no obvious correlated variability can be seen between the 37 GHz flares and the other studied bands. There are a few possible explanations for this.

First, the cadence of the Swift observations might be too low to catch the trends seen in the 37 GHz light curve. Assuming that the flare is quasi-simultaneous in different bands, and given the extremely fast timescales in radio, the X-ray observations would likely need to be simultaneous with the radio to catch the X-ray flare. On the other hand, in the case of SDSS J164100.10+345452.7, the X-ray emission stayed elevated for 10-20 days after the radio flare [70]. This might point to a different origin or possibly environment between the two sources. Another option is that the flares at shorter wavelengths are too short to catch with the long Swift exposures, yielding averages with quite similar, low cps values. Considering the faintness of the source, higher-sensitivity observations would be required to achieve higher temporal resolution. Third, it is also possible that there is a delay between the short-wavelength and radio emission, and the high-activity period in the optical-to-X-ray bands did not occur during our observations. Lastly, the radio flares may not always, or in this source, have multiwavelength counterparts. In blazars, orphan flares are seen over the electromagnetic spectrum, including radio frequencies [31, 86]. However, it is hard to say whether the comparison with J1228+5017 makes sense when the source does not seem to host standard relativistic jets.

Nature of the flares

The derived timescales are upper limits, but the flares cannot be arbitrarily short and pulse-like either; otherwise, they would not result in detections at MRO, and these pulses would also be seen in the temporally resolved MRO data. Unfortunately, due to sensitivity limitations, the MRO observations cannot be split into shorter integrations to study possible flux density changes within one $\sim1800\ \mathrm{s}$ integration.

The most extraordinary — and constraining – properties of these flares are their extremely short rise and decay timescales coupled with very high variability amplitude. Such timescales have not been seen in AGN before; excluding sources likely affected by interstellar scintillation, the blazar S5 0716+714 might come closest. Based on intraday variability, the $T_{\mathrm{b,var}}$ of S5 0716+714 has been estimated to possibly be as high as $\sim10^{17}\ \mathrm{K}$ (e.g., [84]), but the inferred rise and decay timescales are still of the order of days, and thus much longer than what is seen in J1228+5017 [61, 67]. With this combination of properties, these extreme flares do not seem to have similar counterparts in the literature, and they might represent a new kind of AGN variability.

The possibility that these flares could be caused by interstellar scintillation was rejected by Järvelä et al. [35], but is worth a brief discussion here. The Galactic coordinates of J1228+5017 are $l = 132.0$, $b = +66.4$, so the column density of free electrons can be expected to be very low [26]. The observing frequency of 37 GHz is also well in the weak-scattering regime, which usually produces amplitude variability of only 1–2%. Moreover, to undergo scattering, the angular diameter of the source has to be approximately the same size or smaller than the first Fresnel zone of the scattering screen [47]. At 37 GHz, this is $\sim10\ \mu\mathrm{as}$; based on VLBA observations of these sources, no mas-scale components were detected.

Origin of the flares

With the currently available data, we cannot conclusively determine the physical origin of the radio flares. Revisiting the three hypotheses presented by Järvelä et al. [35], based on the timescales derived here, we can safely rule out the jet/jet base interaction with a BLR cloud/star, which is predicted to produce flares lasting 1-2 days. Next, we will discuss one of these hypotheses — magnetic reconnection — as a possible physical origin of these flares.

One viable explanation is magnetic eruptions near the black hole. A near-Eddington accretion flow carries magnetic loops of alternating polarity onto the horizon, where each loop is annihilated by the next one of opposite sign [64, 14]. With $r_{\mathrm{g}} = GM/c^2 \approx1.3 \times10^{12}\ \mathrm{cm}$ and $\dot{M} \approx L_{\mathrm{Edd}}/0.1c^2 \approx10^{25}\ \mathrm{g\ s^{-1}}$, the horizon mass density $\rho\approx\dot{M}/4\pi r_{\mathrm{g}}^3 c \approx2 \times10^{-11}\ \mathrm{g\ cm^{-3}}$ lets a loop saturate at an equipartition magnetic field strength of $B_{\mathrm{H}} \approx(8\pi\rho c^2)^{1/2} \approx6 \times10^5\ \mathrm{G}$, so a loop of size $h$ stores

A burst radiates an isotropic-equivalent $E_{\mathrm{rad}} \sim10^{45}$–$10^{47}\ \mathrm{erg}$ per hour ($\nu L_{\nu} \sim4.5 \times10^{43}\ \mathrm{erg\ s^{-1}}$ at the 0.6 Jy peak), less than what one loop stores even without beaming. The loops also arrive on the inflow time of the magnetically elevated inner disk, $t_{\mathrm{in}} \approx$ days (e.g., [5]), matching the days-long intervals between bursts.

Reconnection between the loops can eject the annihilated flux into the polar cavity. The cavity contains mostly diluted $e^{\pm}$ plasma created by a pair cascade near the corona (e.g., [54]), so the region has a high magnetization $\sigma= B^2/4\pi n_{\pm}m_{\mathrm{e}}c^2 \gg1$ and accelerates impulsively to $\Gamma\gg1$. An observer on the axis sees $\delta\approx2\Gamma$; $\delta\approx140$ brings the intrinsic brightness temperature of the flare, $T_{\mathrm{b,var}}/\delta^3 \approx10^{12}\ \mathrm{K}$ ($3 \times10^{10}\ \mathrm{K}$ for the 8 h burst duration), to the inverse-Compton limit [37]. Induced Compton scattering absorbs radiation of this brightness temperature in any plasma it crosses: the corona and the disk wind are opaque to the radio emission, and the cavity pairs become transparent only far away, beyond $r \gtrsim10^3\ r_{\mathrm{g}} \approx10^{15}\ \mathrm{cm}$, so the shell must emit the burst beyond this radius and it can be viewed along the axis through a channel free of broad-line clouds.

The shell can dissipate in several ways. Internal reconnection and turbulence between its stripes of alternating polarity energise the pairs to $\gamma\sim\sigma$; as they cool, pairs with $\gamma\sim750$ ($\gamma\approx15$ for the 8 h burst) would emit at the 37 GHz band at $r \sim4 \times10^7\ r_{\mathrm{g}}$ ($1.6 \times10^{14}\ r_{\mathrm{g}}$), where the comoving field is $B' \approx2 \times10^{-4}\ \mathrm{G}$ (0.5 G). This synchrotron emission is broadband and linearly polarised up to $\sim70\%$. It could be accompanied by a simultaneous hard X-ray flare from the energised pairs (e.g., [57]).

Alternatively, the shell keeps accelerating as $\Gamma\propto r^{1/3}$ to $\Gamma\sim10^3$ and drives a magnetised shock into the ion-loaded disk wind at $r \gtrsim10^4\ r_{\mathrm{g}}$, where the reduced wind magnetisation $\sigma_{\mathrm{w}} \sim0.2$ lets the synchrotron maser radiate the shock precursor at an efficiency of $10^{-3}$ [74]. This emission is coherent, linearly polarised (in the X mode), and narrow-band ($\Delta\nu/\nu\approx0.6$), which reproduces the 15 GHz non-detection, and drifts downwards in frequency as the shell decelerates in the wind. Simultaneous multi-frequency and polarimetric radio coverage of one burst, together with a hard X-ray ($>50\ \mathrm{keV}$) exposure, could distinguish the two channels.

Another observable expected to be associated with magnetic reconnection events in AGN coronae is the Neupert effect (e.g., [56, 24, 46]), whereby the light curves in both hard X-rays ($L_{\mathrm{HX}}(t)$) and radio ($L_{\mathrm{R}}(t)$) are correlated with the time derivative of the light curve in soft X-rays ($L_{\mathrm{SX}}(t)$),

as observed from the solar corona during some solar flares. In the solar scenario, a stream of charged particles is accelerated along a closed magnetic field loop and rejoins the chromosphere at the loop footprint (e.g., [28]). Some of these particles may escape the flow as a coronal mass ejection (CME). The hard X-ray and radio flares are attributed to the instantaneous nonthermal bremsstrahlung and synchrotron emission of the downward-streaming particles. The more slowly variable soft X-rays are attributed primarily to thermal bremsstrahlung as the particle stream deposits the bulk of its energy into the foot of the loop, heating the chromospheric plasma and evaporating it into the corona. Vardi et al. [80] report the observation of a solar-like AGN flare event in NGC 3783, wherein the hard and soft X-ray variability appears consistent with the Neupert effect from a magnetic reconnection event, and they suggest that an ultra-fast outflow observed by Gu et al. (2025) during the same flare event is a CME-like ejection from the magnetic loop stream. They also find a trailing rise in the UV, which they attribute to hard X-ray photons from the flare illuminating the accretion disk and enhancing its UV emission.

With the present data, we cannot confirm or reject such an event for J1228+5017; while we do not see corresponding X-ray or UV enhancement in the light curves (Figure 2), it may simply be that we do not have the necessary observing cadence. On the other hand, the Güdel-Benz relation should also hold during CMEs, suggesting that a $10^{4}$ Jy X-ray flare should accompany a 1 Jy radio flare; this kind of X-ray flare would surely have been picked up by X-ray monitoring instruments, such as the Swift/Burst Alert Telescope (BAT, [4]) or the Monitor of All-sky X-ray Image (MAXI, [52]). Thus, it does not seem likely that these flares could be attributed to CME-like phenomena.

Conclusions

In this paper, we investigated the multiwavelength behaviour of the NLS1 galaxy J1228+5017 during a period of high activity at 37 GHz. We combined data from MRO, OVRO, JVLA, and Swift/XRT and UVOT to further constrain the observational properties and the physical origin of these extreme high-frequency radio flares. Our conclusions are the following:

i) The extremely short inferred timescale (e-folding time of 1.1 h) coupled with the high, Jy-level amplitude variability and the absence of persistent relativistic jets likely rules out known AGN radio variability mechanisms, strongly suggesting that alternative explanations are required.

ii) The lack of multiwavelength trends might indicate that these flares are exclusively detectable in the radio regime. Alternatively, their higher-energy counterparts might occur at bands not probed here, such as the hard X-ray band, as suggested by the magnetic reconnection scenario.

iii) The emitting region size, a few mpc, seems to be inconsistent with the expected size of the radio corona at GHz frequencies. Interestingly, the sizes derived from the variability timescale and the low-state flux density are similar, possibly indicating that low- and high-state emission originates from the same region.

iv) Radio emission powered by magnetic eruptions stays a viable hypothesis for the physical origin of the flares. This would indicate high magnetic activity in these sources, as predicted by simulations of near-Eddington accreting supermassive black holes.

While these new data further constrain the properties of the extreme radio flares, their physical origin remains unclear. Targeted observations will be needed to test, for example, the predictions of the magnetic reconnection hypothesis. Future observations should focus on measuring the polarisation of the radio emission as well as obtaining high-time-resolution broadband dynamic radio spectra. It could provide a direct test between broadband synchrotron variability and the narrow-band, frequency-drifting emission expected in the synchrotron-maser scenario. These kinds of observations would tremendously help us narrow down the options. Furthermore, simultaneous radio and hard X-ray observations would directly test the magnetic eruption scenario. Due to the sporadic nature of these flares and their short timescales, these observations will not be easy, but considering that we might be witnessing a totally new type of AGN variability, they will be worth it.

Acknowledgements

We want to thank an MRO observer, I. Syrjärinne, for obtaining some of the MRO data presented here.

JN acknowledges funding by the Research Council of Finland grant No. 372879 (OVERFLOW).

This publication makes use of data obtained at Metsähovi Radio Observatory, operated by Aalto University in Finland.

The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. CIRADA is funded by a grant from the Canada Foundation for Innovation 2017 Innovation Fund (Project 35999), as well as by the Provinces of Ontario, British Columbia, Alberta, Manitoba and Quebec.

This research has made use of data from the OVRO 40-m monitoring program ([69]), supported by private funding from the California Institute of Technology and the Max Planck Institute for Radio Astronomy, and by NASA grants NNX08AW31G, NNX11A043G, and NNX14AQ89G and NSF grants AST-0808050 and AST-1109911.

This research has made use of the NASA/IPAC Extragalactic Database (NED), which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration. This research has made use of the SIMBAD database, operated at CDS, Strasbourg, France.

No AI tools were used in the creation of this manuscript.

Data Availability

The JVLA (Legacy ID: AJ442) data are publicly available in the NRAO Data Archive: https://data.nrao.edu. The MRO data will be made available via CDS. The OVRO 40m light curve underlying this article will be shared on reasonable request to the OVRO 40m Monitoring collaboration. The data underlying this article are publicly available from the Swift archive.

Observation details

The details of the Swift/XRT observations are given in Table A.1.

Sequence	MJD	Start time (UT) (yyyy-mm-dd hh:mm:ss)	End time (UT) (yyyy-mm-dd hh:mm:ss)	Exposure (s)
00016062001	60094.76896	2023-05-30 15:09:17	2023-05-30 21:45:17	5661
00016062002	60097.67922	2023-06-02 05:07:15	2023-06-03 03:28:53	4418
00016062003	60098.51066	2023-06-03 08:05:48	2023-06-03 16:24:53	4403
00016062004	60103.59945	2023-06-08 13:30:31	2023-06-08 15:15:53	2146
00016062005	60104.29692	2023-06-09 03:49:18	2023-06-09 10:22:53	1647
00016062006	60106.25109	2023-06-11 05:01:17	2023-06-11 07:01:52	2119
00016062007	60108.36176	2023-06-13 07:47:57	2023-06-13 09:33:54	2131
00016062008	60110.28413	2023-06-15 05:52:22	2023-06-15 07:45:54	2397
00016062009	60113.83124	2023-06-18 19:53:03	2023-06-18 20:00:55	471
00016062010	60120.64475	2023-06-25 12:11:59	2023-06-25 18:44:52	1918
00016062011	60123.61480	2023-06-28 12:55:43	2023-06-28 16:34:53	2480
00016062012	60126.45835	2023-07-01 09:15:10	2023-07-01 12:44:52	2517
00016062013	60129.53083	2023-07-04 05:24:28	2023-07-04 20:04:19	2379