Unprecedented extreme high-frequency radio variability in early-stage active galactic nuclei

We report on the discovery of one of the most extreme cases of high-frequency radio variability ever measured in active galactic nuclei (AGN), observed on timescales of days and exhibiting variability amplitudes of three to four orders of magnitude. These sources, all radio-weak narrow-line Seyfert 1 (NLS1) galaxies, were discovered some years ago at Aalto University Mets\"ahovi Radio Observatory (MRO) based on recurring flaring at 37 GHz, strongly indicating the presence of relativistic jets. In subsequent observations with the Karl G. Jansky Very Large Array (JVLA) at 1.6, 5.2, and 9.0~GHz no signs of jets were seen. To determine the cause of their extraordinary behaviour, we observed them with the JVLA at 10, 15, 22, 33, and 45 GHz, and with the Very Long Baseline Array (VLBA) at 15 GHz. These observations were complemented with single-dish monitoring at 37 GHz at MRO, and at 15 GHz at Owens Valley Radio Observatory (OVRO). Intriguingly, all but one source either have a steep radio spectrum up to 45 GHz, or were not detected at all. Based on the 37 GHz data the timescales of the radio flares are a few days, and the derived variability brightness temperatures and variability Doppler factors comparable to those seen in blazars. We discuss alternative explanations for their extreme behaviour, but so far no definite conclusions can be made. These sources exhibit radio variability at a level rarely, if ever, seen in AGN. They might represent a new type of jetted AGN, or a new variability phenomenon, and thus deserve our continued attention.


INTRODUCTION
Approximately 10 per cent of active galactic nuclei (AGN) are capable of launching and maintaining relativistic jets (Padovani 2017).Traditionally, these jetted AGN have been often identified using the radio loudness parameter 1 as a proxy for the jet activity: all the jetted AGN were believed to be found among the radio-loud population.Whereas the radio loudness parameter might still serve a purpose when considering bright, high-redshift AGN with steady, powerful jets, and negligible host galaxy contribution, recent studies have ★ E-mail: ejarvela@ou.edu† Dodge Family Prize Fellow in The University of Oklahoma 1 Radio loudness parameter, R, is defined as the ratio between 5 GHz flux density and optical  band flux density.Sources with R > 10 are considered radio-loud, R < 10 radio-quiet (Kellermann et al. 1989).
shown that it utterly fails when faced with the true diversity of AGN jet phenomenon and variability (Padovani 2017;Lähteenmäki et al. 2018).This is especially problematic in the local Universe, where we are able to detect also lower-power jets and outflows in AGN, and where the host galaxy can have a major contribution to the lowfrequency radio emission, such that disentangling different sources of radio emission poses a problem (Caccianiga et al. 2015;Järvelä et al. 2017Järvelä et al. , 2022)).This can lead to AGN with low-power relativistic jets to be classified as radio-quiet, or non-jetted AGN with strong star formation to be classified as radio-loud (Caccianiga et al. 2015), making radio loudness a problematic proxy for the jet power and activity.
Especially one class of AGN, the narrow-line Seyfert 1 (NLS1) galaxies, have played a major role in revealing the diversity seen in AGN activity, and have revolutionised some long-standing assumptions held about AGN.NLS1s are identified based on the optical spectrum: the full-width at half maximum (FWHM) of their broad H emission line is <2000 km s −1 , and their [O III] emission is weak compared to the broad H: ([O III])/(H)<3 (Osterbrock & Pogge 1985;Goodrich 1989).They often also exhibit strong Fe II emission, confirming the unobstructed view of the central engine.The narrow FWHM(H) can be attributed to low rotational velocity around a low-mass supermassive black hole (10 6 -10 8  ⊙ , Peterson 2011; Komossa et al. 2018).The low-mass hypothesis is supported by reverberation mapping studies (Wang et al. 2016;Du et al. 2018), predominantly turbulence-dominated Lorentzian emission line profiles (e.g., Kollatschny & Zetzl 2011;Sulentic et al. 2000;Berton et al. 2020a), the existence of tidal disruption events in NLS1s (e.g., Frederick et al. 2021), and the prevalence of disk-like host galaxies with pseudo-bulges (e.g., Järvelä et al. 2017;Olguín-Iglesias et al. 2020;Varglund et al. 2022).The luminosities of NLS1s, comparable to those of higher black hole mass AGN, such as broad-line Seyfert 1 (BLS1) galaxies, combined with their lower black hole masses indicate that a considerable fraction of NLS1s are accreting close to or even above the Eddington limit (Boroson & Green 1992).This ensemble of properties has led to the conclusion that they are fastgrowing, early-stage AGN (Mathur 2000), possibly experiencing one of their first activity cycles.
Based on their properties, NLS1s were not expected to show prominent jet activity, as the ability to launch and maintain powerful relativistic jets was considered to be exclusively a property of massive elliptical galaxies, hosting the most massive black holes (Laor 2000).However, contradictory to this jet paradigm several NLS1s were found to exhibit blazar-like properties in radio band (Komossa et al. 2006;Yuan et al. 2008), and finally the first NLS1 was detected at gamma-rays -indisputably produced by relativistic jets -in 2009 (Abdo et al. 2009).Since then ∼20 NLS1s have been detected at gamma-rays (Romano et al. 2018;Paliya 2019), and several dozen new candidates have been identified (Foschini et al. 2021(Foschini et al. , 2022)).Furthermore, additional ∼50 NLS1s have been confirmed to host jets via radio imaging (e.g., Richards & Lister 2015;Lister et al. 2016;Berton et al. 2018;Chen et al. 2020Chen et al. , 2022)).NLS1s with relativistic jets share similar properties with the non-jetted NLS1 population, and thus broke the jet paradigm beyond any doubt.These jetted NLS1s are also the first AGN with systematically high Eddington ratios to host relativistic jets.Blazars, in general, have Eddington ratios < 0.1 (Heckman & Best 2014), and it was believed that AGN with Eddington ratios significantly higher than that are very rarely capable of launching jets, though some exceptions exist (e.g., Belladitta et al. 2022).Recently, also general relativistic (radiative) magnetohydrodynamics (GRRMHD) simulations have shown that efficient and powerful collimated jets are formed in systems with high Eddington ratios, even exceeding unity, if the state of magnetically-arrested accretion (MAD) is reached (McKinney et al. 2017;Liska et al. 2022).Thus it seems that our earlier beliefs regarding relativistic jets were mainly a product of observational biases, for example, concentrating the studies only on the brightest or radio-loudest AGN.It has been suggested that jetted NLS1s represent an early stage of the evolution of jetted AGN, and that they will eventually grow into flat-spectrum radio quasars (FSRQ) and radio galaxies (Foschini et al. 2015;Berton et al. 2017).If this is the case, they offer us an unprecedented opportunity to study the very first stages in the evolution of powerful AGN with relativistic jets.
Intriguingly, the radio properties of NLS1s are very diverse: only 15 per cent of them have been detected in radio (Komossa et al. 2006;Järvelä et al. 2015), and include a continuum of sources from host-dominated to relativistic jet-dominated (Järvelä et al. 2022), whereas the majority of 85 per cent seem to be totally radio-silent.
However, NLS1 samples often suffer from misclassifications, and include a significant fraction of BLS1s and intermediate-type AGN that affect the population-wise statistics.Indeed, an ongoing investigation utilising a carefully selected sample of NLS1s and new radio surveys, such as the LOw-Frequency ARray (LOFAR) Two-metre Sky Survey (LoTSS) and the National Radio Astronomy Observatory (NRAO) Very Large Array Sky Survey (VLASS), indicates that the radio detection fraction among NLS1s is even lower, around ∼8 per cent (Varglund et al. in prep.).To understand the nature of this seemingly heterogeneous class and how different NLS1s are related, it is necessary to study the population as a whole.Most studies have concentrated on the most obvious radio-bright NLS1s, whereas the radio-faint and -silent population has been scarcely investigated.

The road so far
A different approach was adopted at the Aalto University Metsähovi Radio Observatory (MRO, Finland), where several hundreds of jetted AGN are frequently monitored at 37 GHz.In addition to the usual suspects, that is, NLS1s that are bright in radio, two samples of NLS1s were selected for monitoring based on totally distinct criteria, independent of their radio properties.One sample consisted of NLS1s residing in very dense large-scale (Mpc-scale) environments, such as superclusters (Järvelä et al. 2017), and the other was compiled from NLS1s exhibiting spectral energy distributions (SED) that seemed favourable for 37 GHz observations.Eight NLS1s from these samples, four from each, were detected at flux density levels of several hundred mJy (Lähteenmäki et al. 2018).What makes these sources extraordinary is that most of them had been deemed to be radio-silent or had only very faint previous radio detections.Seven sources have been detected several times, strongly suggesting that these are genuine detections of recurrent radio flares in these sources.The most likely emission mechanism to produce such high-amplitude, rapid variability at a radio frequency this high is the synchrotron emission of a relativistic jet (however, see Sect. 5).Additional evidence was obtained when one of the sources was identified as a new gamma-ray emitter, and has since been seen brightening in X-rays soon after an MRO-detected flare (Romano et al. 2023).
Only two of these sources had previous radio detections, and only at mJy levels in the Faint Images of the Radio Sky at Twenty-Centimeters survey (FIRST) and the NRAO Very Large Array (VLA) Sky Survey (NVSS), while the rest were non-detections.To decipher this puzzling behaviour and to discriminate between the different hypotheses of their nature, the sources with several MRO detections were observed with the Karl G. Jansky VLA (JVLA) in Aconfiguration in L, C, and X bands.Instead of clarifying the situation, these observations raised more questions.Two of the sources were non-detections and the remaining sources had flux densities ranging from a few tens of Jy to a few mJy, all of them consistently showing steep spectra below 9 GHz (see Fig. 6 in Berton et al. 2020b).Three of them showed slightly extended radio morphology.In a closer inspection, exploiting spatially resolved spectral index maps, it was found that at least two of these sources show signs of flat core spectrum (Järvelä et al. 2021) and thus the presence of a partially optically thick radio core.The JVLA and the MRO observations are not simultaneous, but such an extreme, similar behaviour observed in several sources indicates that it is real, not just a curiosity.
However, the beam size of MRO (∼2 arcmin) is considerably larger than the beam size of the JVLA in A-configuration (∼arcsecscale).It is therefore important to consider the possibility that the discrepancy between the flux densities of the JVLA and MRO could arise from different beam sizes.This seems improbable when taking into account the properties of the emission.Due to the redshift of these sources the JVLA observations probe kpc-scale structures.The angular sizes of these sources in the optical band are between 2 and 12 arcsec, so we were able to see the whole galaxy in the JVLA observations, in which the smallest field of view -at 9 GHz -was 4.7 arcmin.It is hard to explain such strong and variable radio emission in the outskirts of, or even outside, a galaxy.Due to the rapid variability, indicating a small emitting region, it is highly improbable that resolved-out structures could be responsible for this emission.Furthermore, contamination by nearby sources was ruled out in Lähteenmäki et al. (2018).It can thus be assumed that the JVLA and MRO probe the same phenomenon.The effects of different beam sizes is further discussed in Sect.5.2.1.
Since the low frequency flux densities are consistent with FIRST there is no need to assume that these NLS1s have undergone drastic changes, for example, triggering of jets, but it cannot be ruled out either.Thanks to the MRO data we know these sources most likely host relativistic jets, but their radio emission below 9 GHz (X band) seems to be consistent with star formation, with little or no contribution from the AGN.Extrapolating, or even assuming a flat radio spectrum up to 37 GHz would mean that in the quiescent state the flux density would be less than a mJy, which, in the most extreme case, would require a nine-thousand-fold increase during flares.This would be very extreme, and a more plausible explanation is that the spectrum turns inverted at some point above 9 GHz, as indicated by the MRO data.This kind of behaviour is commonly seen in kinematically young AGN, for example, high-frequency peakers and gigahertzpeaked sources (O'Dea & Saikia 2021).In these sources the convex radio spectrum is explained by synchrotron self-absorption (SSA) in a young, parsec-scale jet.However, even in these sources the peak frequency does not usually exceed 10-15 GHz, which in contrast seems to be the case in our sources.
An alternative to SSA could be free-free absorption (FFA), which also allows more inverted spectral indices than SSA (Rodriguez et al. 1993), requiring less extreme variability at 37 GHz.Some cases where the turnover frequency stays consistently high have been found (tens of GHz, Doi et al. 2016), and usually this behaviour is explained by FFA.This could be the case also in these NLS1s: if these sources are kinematically young AGN, FFA could happen in the shocked ionised ambient clouds in front of the jet head (O'Dea & Saikia 2021).Alternatively, the required ionised gas could be provided by the enhanced circumnuclear star formation activity often seen in NLS1s (Sani et al. 2010;Winkel et al. 2022).Either way, these NLS1s with jets that are almost totally absorbed at low radio frequencies are the last nail to the coffin of the radio loudness parameter as a universal proxy for the jet activity of AGN, and urge us to expand our horizons when it comes to our understanding of the diversity of AGN jets.
To discern between these alternatives, we observed seven of these sources with the JVLA in X, Ku, K, Ka, and Q bands.These observations were complemented by Very Long Baseline Array (VLBA) observations at 15 GHz, and single-dish observations at 15 and 37 GHz, using the OVRO 40m telescope and the MRO telescope, respectively.In Sect. 2 we introduce the sample, in Sects.3.1 through 3.4 we describe the performed observations, and the data reduction and analysis, in Sect. 4 we present our results, in Sect. 5 these results and their implications are discussed, and in Sect.6 we provide a brief summary of this work.Throughout this paper, we adopt a standard ΛCDM cosmology, with a Hubble constant  0 = 72 km s −1 Mpc −1 , and Ω Λ = 0.73.

SAMPLE
The sample includes seven radio-weak NLS1s repeatedly detected at Jy-level flux densities at 37 GHz at MRO.The eighth such source was dropped because it was detected only once.Originally these sources were selected for the MRO AGN monitoring based on their dense large-scale environments (Järvelä et al. 2017) or SEDs that suggested that they could be detectable at high radio frequencies (Järvelä et al. 2015).The basic properties of the sample are summarised in Table 1.These sources are very similar to the general NLS1 population: all have a black hole mass less than 10 8  ⊙ (Järvelä et al. 2015;Lähteenmäki et al. 2018), and six of them are hosted in a disk-like host galaxy (Järvelä et al. 2018;Olguín-Iglesias et al. 2020;Varglund et al. 2022), whereas the morphology of the highest- source is unknown.

Observations and pre-processing
We observed our sample with the JVLA in A-configuration in five different bands, X, Ku, K, Ka,and Q,centred at 10,15,22,33,and 45 GHz,PI Järvelä).The dates and integration times of the JVLA observations are given in Table A1.The total bandwidth was 4 GHz in X, 6 GHz in Ku, and 8 GHz in K, Ka, and Q band, each band divided to 128 MHz subbands, consisting of 64 channels of 2 MHz.The NLS1 (Berton et al. 2017) 3C 286 was used as the bandpass and flux density calibrator for each source, and each source had an individual nearby, bright source that was 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 7, 7, 12, 12, and 25 Jy beam −1 in X, Ku, K, Ka, and Q, respectively.We were able to reach these levels in most cases.
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.In addition, the data were checked manually and any remaining bad data were flagged, producing the SRDP measurement set for each source.We also re-checked all the data manually, but no additional flagging was required.In further data processing and analysis we used the Common Astronomy Software Applications (CASA) version 6.2.1-7.We split the data of our sources from the measurement set separately in each band averaging over time (timebin = 10 s) and frequency (width = 64, to average 64 channels to form one output channel per subband).Before the actual imaging of the targets, we produced radio maps of the size of 2.7 arcmin × 2.7 arcmin, or the whole primary beam, to check the whole beam of the MRO and OVRO telescopes to identify any other sources of radio emission within them.We did not find other strong radio emitters in any of these fields, further supporting the assumption that the radio emission detected at MRO is coming from the NLS1 nucleus.

Radio maps and measurements
We used the tclean algorithm with interactive cleaning in CASA to produce the radio images of our sources.The cell size was chosen so that the synthesised beam is properly sampled, meaning a cell size of 250, 150, 100, 70, and 50 mas in X, Ku, K, Ka, and Qbands, respectively.The image size was chosen so that the whole galaxy fits into the image, taking into account the varying cell sizes in different bands, and the redshifts of our sources.We used Briggs Table 1.Basic properties of the sample.Columns: (1) source name in the SDSS, the superscript indicates the band the coordinates are from, G stands for Gaia; (2) short name; (3, 4) right ascension and declination (J2000); (5) redshift; (6) scale at the redshift of the source; (7) logarithmic black hole mass, taken from Lähteenmäki et al. (2018); (8) large-scale environment, taken from Järvelä et al. (2017); (9) host galaxy morphology, PB = pseudo-bulge, taken from  Järvelä  et al. (2018),  Olguín-Iglesias et al. (2020),  Varglund et al. (2022).weighting, with robust = 1.8, in all cases.Some sources appear to be slightly hexagonal (for example, J1522+3439), possibly due to the sidelobes.In these cases we trialled with robust values closer to uniform weighting to suppress the sidelobes but there was no visible difference, so we decided to maximise the sensitivity and use the same robustness value for all sources.No source was bright enough to be self-calibrated.We used the mtmfs deconvolver with nterms = 2 and scales = 0 in case some sources would be bright and extended enough to produce spatially resolved in-band spectral index maps, which turned out not to be the case.However, due to this we did the wide-band primary beam correction separately with widebandpbcor.

SDSS
We fitted each detected source using a 2D Gaussian to obtain the central coordinates and the peak flux density and its error.In cases of extended sources we measured the emission inside 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.In case the source was not detected, we report 3 upper limits.The results are given in Sect.4.1-4.7,and the radio maps shown in App.B.

Observations
We observed our sample also on milliarcsecond scale using the VLBA in the Ku band, centred at 15.1 GHz (Project BJ 109, PI Järvelä).The observations were carried out during one 10 hr long experiment on 2022-02-08.The recording setup used the Digital Downconverter (DDC) system of the Roach Digital Backend (RDBE) with four 128 MHz wide subbands -giving a total bandwidth of 512 MHz -two circular polarisations, and two-bit sampling, resulting in a total recording rate of 4 Gbps.
Due to the potentially low compact flux densities of the target sources, the observations were carried out using the standard phasereferencing technique, i.e. a rapid switching between the target and a nearby calibrator.The phase-reference calibrators together with their distances from the targets, their VLBI scale flux densities, and the used source-switching duty cycles are given in Table A2.Each target source had 38 min total on-source integration time.Bright flat-spectrum radio quasars 3C 279 and 3C 345 were observed for two 5 min long scans and for three 3 min long scans, respectively.They were used as fringe finders and, more importantly, as calibrator sources for determining instrumental delays and bandpass shapes.Nine out of ten VLBA antennas participated in the observations, since Hancock was out due to a frozen focus/rotation mount.

Data reduction
The recorded station data were correlated with the VLBA DiFX correlator in Socorro using 0.5 MHz wide spectral channels and 1 s correlator integration time.This allows a relatively wide field of view, >4" from the phase centre, to be searched for compact sources.
The data were calibrated in the Astronomical Image Processing System (AIPS; Greisen 2003) using standard procedures for phasereferencing observations.The calibration started with a priori corrections to the station parallactic angle, updates to the Earth Orientation Parameters, and first-order removal of dispersive ionospheric delays using total electron content maps derived from the Global Navigation Satellite System (GNSS) data.Instrumental delays and phase offsets between subbands were removed by fringe-fitting a single scan of the bright calibrator 3C 279.A priori amplitude calibration included corrections to sampler threshold levels by using autocorrelations, bandpass calibration using again a scan on 3C 279, and conversion of raw correlator coefficients to Janskys by applying measured system temperatures and gain curves.
The phase reference calibrators as well as the bright calibrators 3C 279 and 3C 345 were fringe-fitted using the AIPS task fring and combining subbands and using an integration time of either 2 min or the scan length, whichever was shorter.The fringe-fitting gave excellent results; the percentage of failed solutions was typically ∼1 per cent.The fringe-fitting solutions from the phase-reference calibrators were applied to both the calibrators and the target sources.The relative R-L delays were corrected by cross-hand fringe-fitting of a single scan of 3C 279.After this step, we imaged the calibrator data in Difmap (Shepherd 1997) and loaded the images back to AIPS.The calibrator images were used to derive phase self-calibration solutions for the calibrator data using the AIPS task calib and 10 s integration times.These phase solutions were then applied to the target sources.As the last correction, we also used the amplitude self-calibration solutions from imaging the bright calibrators 3C 279 and 3C 345 to fine-tune the amplitude calibration for those antennas and subbands that had an average amplitude self-calibration solution deviating more than 5 per cent from unity.After this step, the target data were ready for imaging.

Imaging and searching for the target sources
While we had quite accurate a priori positions of the target sources based on the previous JVLA data (positional uncertainties less than 10 mas), we still wanted to search for an area that covers most of the galaxy in case the variable emission seen in the single-dish data does not come from the JVLA core.To achieve this, for each target source we generated a set of naturally weighted images with a field-of-view of 820×820 mas that covered an area of 7.4"×7.4"centred on the JVLA position using the multifield option of the AIPS task imagr.
The image rms was ∼ 60 Jy/beam for all the target sources which is at the expected thermal noise level.Since we searched for a large area covering one million synthesised beam areas per image, we set the detection threshold to 6 to avoid picking noise spikes.No sources were detected, and in Tables 3, 5 , 8, 11, 13, 15, and 18 we quote 6 upper limits for the VLBA data.

Single-dish data
In addition to radio interferometric data, we obtained nonsimultaneous single-dish monitoring data for all of these sources from MRO and OVRO; these data will be published here.We also have 1-3 epochs of single-dish observations per source from the Effelsberg 100-m Radio Telescope between 4.5 and 45 GHz, and one epoch of 2 and 1.15 mm observations with the New IRAM Kids Arrays (NIKA2) instrument on the Institut de Radioastronomie Millimétrique (IRAM) 30-m Radio Telescope on Pico Veleta for five sources.The Effelsberg and IRAM data, complemented by MRO and OVRO data from the same time period, will be published in an upcoming paper.

Metsähovi Radio Observatory
The measurements included in this study are part of the large ongoing AGN monitoring programme at 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 used 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 a S/N < 4 are handled as non-detections.(2011) and more details specific to the NLS1 observations are given in App.A4.These seven sources were added to the OVRO AGN monitoring programme in July 2020, and since then three of them have been detected with S/N > 4.This paper includes OVRO data until June 2022.

Archival data
In addition to the new data obtained we also used already published data of these sources.We included the JVLA A-configuration L, C and X band data from Berton et al. (2020b) taken in September 2019.We also included the LOFAR LoTSS Data Release 2 (DR2) data with a central frequency of 144 MHz (band 120-168 MHz) (Shimwell et al. 2022).All of our sources reside within the published region of the sky.The resolution of LOFAR LoTSS DR2 is 6 arcsec, the median rms sensitivity is 83 Jy beam −1 , the flux density scale accuracy is ∼10 per cent, and the astrometric accuracy is 0.2 arcsec.We used a 1.2 arcmin search radius to check the whole MRO beam area.In addition, we checked the Stokes I continuum radio maps to correctly identify the NLS1, and any other possible radio sources, and to visually cross-match the radio sources with any optical sources.Last, we included NRAO VLASS Epoch 1 and 2 data.The angular resolution of VLASS is ∼2.5 arcsec, and it covers the entire sky north of  = −40 deg.In this paper we use data based on the Quick Look and single epoch imaging, which have a systematic ∼15 per cent underestimation of the flux density values at  peak > 3 mJy beam −1 .We used the same search radius as for LOFAR.These data are discussed in detail in the individual source sections.

RESULTS
The results for each source are given in the following Sections.In addition to the radio map measurements, we calculated the redshiftand k-corrected radio luminosities as: where  is the central frequency of the band in Hz,   the observed flux density in erg s −1 cm −2 Hz −1 ,  2  the luminosity distance in cm, and  the spectral index of the emission.For simplicity we used  = 0 in all calculations.Even drastic changes in  do not significantly affect the luminosity, i.e. the order of magnitude remains the same.Furthermore, since our sources are variable they do not have a characteristic spectral index.The luminosities are given in the Tables in the following Sections for individual sources.
Since our sources are only marginally extended or point-like and their JVLA spectra show a consistent slope throughout the detected bands, it is unlikely that in-band spectral index maps could yield significant new information regarding their spectral properties.Thus we calculated only the traditional spectral indices between new detections with interferometric arrays using both the peak flux densities and the integrated flux densities.
Additionally, we used temporally close consecutive 37 GHz detections to estimate the properties of the flares.The details of the calculations and the results are given in Sect.4.8, but referred to in the following sections for the individual sources.

SDSS J102906.69+555625.2
So far J1029+5556 has been detected at 37 GHz at MRO and at 15 GHz at OVRO (Table 4).It has not been detected in any radio interferometric observations (see Table 3).J1029+5556 has the highest redshift,  = 0.451, in this sample, and due to this it is also the only source that is missing the host galaxy morphology information.Interestingly, it was detected at MRO only three times in 2016-2017, with moderate flux densities around 500 mJy and below, and has not been detected after that, though it has not been observed very frequently in the past few years.Its overall detection percentage at 37 GHz is 6.1 per cent, and the mean luminosity   = 9.5 × 10 43 erg s −1 .The lack of recent detections might indicate a change in the activity level of the nucleus, though it was detected by OVRO in 2020, indicating that the activity has not totally halted.Whether the amplitude of the variability has changed or if the most drastic variability has moved to lower frequencies cannot be determined based on these data.J1029+5556 is not present in LOFAR maps, but there is one radio source, which lacks an optical counterpart, in the LOFAR map within the MRO beam.However, the source is faint, with a flux density of ∼1 mJy, and we do not see signs of it in the JVLA data.The non-simultaneous radio spectrum of J1029+5556 is shown in Fig. 1 and the light curves in Figs.C1 and C8.

SDSS J122844.81+501751.2
J1228+5017 is detected with the JVLA in all other bands except the Q band, and it is also detected by LOFAR at 144 MHz (Table 5).It is not properly resolved in any JVLA band (Figs.B1-B5).In the 144 MHz radio map it seems to be extended toward north-west, but upon closer inspection the extended part turns out to be a nearby galaxy that can also be seen in optical images.The radio spectrum, shown in Fig. 2, has a constant slope of −0.7 from 144 MHz to X band, above which the slope flattens considerably (Table 6).The low-frequency spectral index is consistent with the characteristic star formation activity spectral index of −0.7 and the flux density levels could be explained by star formation (Berton et al. 2020b).Though it should be noted that also the spectral index of optically thin synchrotron emission by shock-accelerated electrons in jets is around −0.7.The spectrum shows the characteristic spectral turnover, or spectral index flattening, toward lower frequencies where the emitting medium starts to become opaque to radio emission (Condon 1992).In principle, the high-frequency spectral index is very close to the thermal free-free emission spectral index of −0.1, that in star forming galaxies has an increasing contribution toward higher frequencies, whereas the steep synchrotron emission from supernovae becomes less important.However, the change in the slope between the nonthermal and thermal emission dominated spectral regions should not be this drastic (Klein et al. 2018).Instead, the flattening could be due to a third component, the flat radio core of the AGN that becomes detectable when the emission produced by star formation weakens.Spatially resolved spectral index maps in L, C, and X bands support this scenario since despite the overall steep spectral index, the core spectral index in these bands is significantly flatter (Järvelä et al. 2021).X band also shows a peak flux density decrease from 0.184 ± 0.008 mJy beam −1 in Berton et al. (2020b) to 0.128 ± 0.005 mJy beam −1 in these observations.The JVLA configuration and the rms of the maps are the same for both observations, but the central frequencies are slightly different (9 vs 10 GHz), thus the difference could be due to the slightly different beam sizes, since the source is partially resolved.
J1228+5017 has been detected at MRO seven times, with the last detection in 2019, and has a detection percentage of 15.2 per cent and a mean luminosity of   = 2.6 × 10 43 erg s −1 .The single-dish detections are listed in Table 7 and the light curves are shown in Figs.C2 and C9.However, the source does not seem to have totally gone into slumber as it has been detected again recently (Järvelä et al. in prep.).

SDSS J123220.11+495721.8
In the earlier JVLA observations J1232+4957 was detected in L and C bands, but not in X band.In the new observations it is also detected in X and Ku bands, but only at a 3 level (Table 8).Also LOFAR detected J1232+4957 at 144 MHz.It remains unresolved in all interferometric observations (Figs. B6 and B7).Its radio spectrum, in Fig. 3, clearly shows a steepening slope toward higher frequencies.The spectral index between 144 MHz and X band is −0.56 ± 0.08, and between X and Ku bands −1.49 ± 0.59.(Table 9).The interferometric flux densities and the spectral properties of J1232+4957 can be explained by star formation activities, and AGN contribution does not seem necessary.
On the other hand, J1232+4957 has been detected at MRO several times with an overall detection percentage of 10.6 per cent.The mean luminosity of the detections is   = 2.8 × 10 43 erg s −1 .The last detection, however, is from 2019 (Table 10).The 37 GHz flux

SDSS J150916.18+613716.7
J1509+6137 is an intriguing source as it has clearly the highest detection percentage at 37 GHz -25.3 per cent -but it has not been detected in any JVLA band.The MRO detections have an average luminosity of   = 2.5 × 10 43 erg s −1 .The light curves are shown in Figs.C4 and C11, and the radio data are given in Tables 11  and 12.The brightest MRO flares exceed 1 Jy, indicating extreme variability of four orders of magnitude.J1509+6137 also has several double-detections within a week from each other.These detection pairs were used to estimate the flare characteristics (Table 21) and are discussed in Sect.4.8.J1509+6137 was not detected by LOFAR, but based on the LOFAR LoTSS DR2 there are two other radio sources within the MRO beam.Neither of these sources have optical counterparts, and both of them are faint, around 0.4 and 0.8 mJy.They are not seen in the JVLA data.J1509+6137 seems to be totally absent in radio -except during the 37 GHz flares -and does not even show detectable amounts of radio emission from star formation.

SDSS J151020.06+554722.0
J1510+5547 has a high detection percentage of 17.6 per cent at 37 GHz (Table 14 and Figs. C5 and C12).It was last detected in 2019 even if the number of annual observations has stayed roughly the same.The mean luminosity at 37 GHz is   = 8.7 × 10 42 erg s −1 .It was detected in L, C, and X bands in our previous JVLA observations, but remained a non-detection in all bands, X through Q, in the recent observations (Table 13).The radio spectrum of J1510+5547 is shown in Fig. 5.The X band upper limit is very close to the earlier X band detection flux density, and considering that the central frequencies of the two observations differ by 1 GHz, it is likely that the recent non-detection is due to the source being very close to the detection limit.
This source is also detected by LOFAR and seems to be marginally resolved.There is another radio source north-east of it and within the MRO beam.This source is faint, has no optical counterpart, and is not seen in any JVLA band.The projected distance between J1510+5547  and the source is more than 40 kpc, thus it is unlikely that it is related to our source.The radio spectrum below 10 GHz is consistent with that of star forming galaxies with the characteristic spectral turnover seen toward lower frequencies.

SDSS J152205.41+393441.3
J1522+3934 is a nearby source ( = 0.077) that resides in a disk galaxy that is merging with a non-active galaxy (Järvelä et al. 2018).It shows almost symmetrical resolved emission on west/north-west and east/south-east sides of the nucleus from 144 MHz to Ku band, and is detected up to Ka band (Table 15 and Figs.B8-B11).Interestingly, the extended radio emission is perpendicular to the host galaxy, indicating that it does not originate from the star formation activity in the host (Järvelä et al. 2021).To explain the 37 GHz flaring in J1522+3934 the jet emission needs to be relativistically boosted and thus the jet needs to point close to our line of sight.If this is the case, the extended emission would be a relic of past activity -unless the jets are very bent, pointing at us close to the nucleus and turning perpendicular at larger distances.The spatially resolved spectral index map in the L band does show regions of steeper spectral index around −1.0, possibly indicative of synchrotron cooling.
The radio spectrum of J1522+3934, in Fig. 6, has a very sta- ble slope around −0.7 from 144 MHz all the way to the Ka band (Table 16).The VLASS points seem to deviate from this which is surprising considering that the Quick Look flux densities should underestimate the real flux densities.Overall the spectrum seems to be consistent with optically thin radio emission and we can assume its predominant origin to be the AGN.J1522+3934 has the record 37 GHz flux density among our sources at 1430 mJy, whereas the other detections are much more modest.Its detection percentage at MRO is only 3.9 per cent, and the mean luminosity is   = 2.7 × 10 42 erg s −1 .In addition to these detections, it has also been detected at 15 GHz at OVRO on three different dates (Table 17), with a maximum flux density of 45 mJy.The light curves of J1522+3934 are shown in Figs.C6 and C13.

SDSS J164100.10+345452.7
J1641+3454 is the only one of our sources with a statistically significant gamma-ray detection (Lähteenmäki et al. 2018), usually considered as proof of the presence of relativistic jets.Interestingly, its detection rate at 37 GHz is the lowest in the sample at 1.5 per cent.Its 37 GHz flux densities are modest, generally around 500 mJy and below, indicating that most of its flaring activity might not exceed the MRO detection threshold.Its average 37 GHz luminosity is   = 9.9 × 10 42 erg s −1 .J1641+3454 has also been detected at 15 GHz at OVRO with a flux density of ∼30 mJy (Table 20).
J1641+3454 was a target of an intense 20-month multiwavelength monitoring campaign in radio, optical, ultraviolet and, X-rays (Romano et al. 2023).During the campaign it flared twice at 37 GHz: the first radio flare was followed by brightening in X-rays, whereas the latter flare was not accompanied by any significant changes at other frequencies.Nevertheless, this was the first detection of a 37 GHz radio flare counterpart at another frequency.
J1641+3454 is detected in X, Ku, K, and Ka bands with the JVLA (Table 18).It is resolved in X and Ku bands, with extended emission seen on the north-west and the south-east sides of the nucleus.This emission is seen also at lower frequencies and it appears to be patchy, which points to star formation rather than the AGN as the origin (Berton et al. 2020b).J1641+3454 is also detected at 144 MHz by LOFAR and at 3 GHz in VLASS.At 3 GHz it is not properly resolved but appears elongated in the north-west/south-east direction similarly to the JVLA maps.Interestingly, in the LOFAR map it seems to be elongated toward south-west.This emission has no optical counterpart, but it is clearly outside the host galaxy of J1641+3454, so it remains unclear whether it is related to J1641+3454.
The radio spectrum of J1641+3454, shown in Fig. 7, clearly has a curvature, and it flattens towards lower frequencies and steepens toward higher frequencies, reaching a spectral index around −1.0.No AGN contribution is required to explain the properties of its high-resolution data radio spectrum and no signs of flattening in the spectrum or in the spatially resolved spectral index maps can be seen (Järvelä et al. 2021).
In addition to the 37 GHz detections, J1641+3454 has also been detected once at 15 GHz by OVRO with a flux density of ∼30 mJy.The OVRO detection is quite close to an MRO detection, within 23 days, but unfortunately in the case of these sources we cannot assume that these detections are necessarily from the same event.However, in the case they were, we can derive a quasi-simultaneous spectral index of 2.70 ± 0.63.Since these detections are not strictly simultaneous and we do not know which stage of the flare the detections represent, the spectral index is only a rough estimate.It agrees with SSA within the errors, but might imply that also another source of absorption is required.The light curves of J1641+3454 are shown in Figs.C7 and  C14.

Flare characteristics using MRO data
We can use the consecutive MRO detections to infer some properties of the radio emission in our sources.Following Valtaoja et al. (1999) and Hovatta et al. (2009) we can estimate the flare rise and decay -folding timescales, variability brightness temperatures, and variability Doppler factors.We performed these calculations for all consecutive detections -that is, there are no non-detections between them -that were less than seven days apart and had different flux densities even when taking the errors into account.We cannot be sure if the two detections are from the same flare, but in case they are not, it means that the variability is even faster and more extreme.We also assume that the maximum amplitude of the flare is equal to the higher of the two flux densities.In case it is not, and the real amplitude of the flare is larger, the timescales would be shorter.Thus, these timescale estimates, and the parameters derived from them can be consider as lower limits.For simplicity, since our knowledge of these sources is so limited, we used the same equation for both rising and decaying flares: where Δ max is the maximum amplitude of the flare in Jy, after subtracting the baseline flux density level,  b ,  max is the epoch of the peak of the flare, and  is the rise or decay time of the flare expressed in days (-folding timescale).We do not know the exact quiescent flux density level, but based on the OVRO observations it cannot be much higher than ∼10 mJy (see Sect. 5.2.1), so we chose this number as the baseline flux density level.The results are shown in Table 21.
In order to estimate the variability Doppler factors of our sources, we calculated the variability brightness temperature,  b, var , (in the where  is the observed frequency in GHz,   is the luminosity distance in metres, and Δ max and  are defined in Eq. 2. The numerical factor corresponds to using  0 = 72 km s −1 Mpc −1 , and Ω Λ = 0.73, and to assuming that the source is a homogeneous sphere.Since estimating the brightness temperature from the flux density variability is based on a causality argument, these values are in fact lower limits.We calculated the variability brightness temperatures for all flares with  values.It should be kept in mind that the brightness temperatures derived from variability are systematically larger by a factor of  2 , where  is the Doppler factor, than those obtained directly from VLBI measurements due to the different dependence on the Doppler factor.
Once we know the variability brightness temperature we can use it to estimate the variability Doppler factor, assuming we know the intrinsic brightness temperature,  b, int : For the intrinsic brightness temperature, we use 5 × 10 10 K (Readhead 1994; Lähteenmäki et al. 1999), which assumes equipartition between the energy densities of the magnetic field and the radiating particles.However, we do not know if these sources really are in equipartition and therefore cannot say how accurate the Doppler factor estimates are.Indeed, the rapid variability suggests that this may not be the case, thus these estimates should be taken with a grain of salt.
Keeping these caveats in mind, the results are reported in Table 21.There are three sources with consecutive MRO detections within one week: J1228+5017, J1509+6137, and 1510+5547, but after excluding all the detections that can be the same within the error bars, only one source, J1509+6137, remains.It has shown two rising and one decaying flare that meet our criteria.In all cases the e-folding timescales are of the order of days, or a maximum of a few weeks, the variability brightness temperatures around 10 14 -10 15 K, and the variability Doppler factors between 5 and 50.These parameters, except the timescale, are comparable to what is seen in flat-spectrum radio quasars (Hovatta et al. 2009).
We can use a simple light travel time argument to infer an approximate size of the radio emitting region.The size needs to be  < /(1 + ).Assuming  ∼ 0, for  of five days this gives 0.0042 pc ×  and for ten days 0.0084 pc × .Taking into account the Doppler factor the size of the emitting region can increase by about an order of magnitude.These sizes are rough estimates since we cannot properly estimate the timescales with the current data, but it is probably safe to assume that the order of magnitude is correct and that the emitting region needs to be milliparsec in size.This indicates that the emission originates close to the black hole, well within the BLR, or from spatially limited regions inside the jet.

DISCUSSION
All of their variability properties considered, these seven sources exhibit flux density variations at a level never observed in AGN before at high radio frequencies.The short variability timescales they show are rare, but not unheard of, even in the radio regime (Rani et al. 2013), whereas the amplitude of the variability -three to four orders of magnitude -coupled with the short timescales, is unprecedented to the best of our knowledge.
Based on the 37 GHz light curves (Figs.C8-C14), including both detections and upper limits (see App. A3 for details), most of the sources are usually detected very close to the detection threshold of MRO.J1509+6137 -that has not been detected in interferometric observations at all -is an exception, and consistently shows activity that is clearly above the detection limit.In general there do not seem to be notable trends in the detections, other than that the sources are detected more when they are observed more, which is not surprising.
In some sources (for example, J1228+5017 and J1232+4957) there seem to be higher upper limits crowding around detections, possibly indicating an increased level of activity during that particular epoch (but see App.A3 for caveats).In others, such as J1641+3454, the detections are embedded amongst upper limits that show no apparent trends of activity.On the other hand, many detections are not accompanied by other nearby observations at all.At OVRO all detections, except the first detection of J1522+3934, are clearly above the detection threshold.However, the detectability at 15 GHz compared to 37 GHz is significantly different.Only three sources have been detected at 15 GHz and the highest detection rate is only 4.5 per cent.The sources with the highest detection rates at 37 GHz have not been detected at 15 GHz at all despite the comparable number of observations.This might indicate that the flaring behaviour is stronger, in terms of the amplitudes, towards higher frequencies.Though, it should be noted that for many sources most MRO detections are from the time before OVRO started monitoring them, so it is also possible that these sources have been less active throughout the OVRO observations.
For some sources (J1509+6137, J1522+3934, and J1641+3454) there are a few MRO detections with OVRO observations within ∼1-5 days before or after the MRO detection.Using these detections and the OVRO upper limits, these quasi-simultaneous observations can be used to estimate a lower limit for the 15-37 GHz spectral index.The spectral index lower limits are around 4 to 5, strongly suggesting that there are external factors in play resulting in the observed phenomenon.
Despite frequent detections at 37 GHz, and some at 15 GHz, all sources were in the low state in the JVLA observations.However, considering the low-to-moderate detectabilities (1.5-25 per cent) and the short timescales of the sources at 37 GHz, it is not infeasible that none of them were flaring at the time of the two epochs of the JVLA observations.
In the following we discuss different phenomena able to cause variability in AGN.It should be kept in mind that the physical explanation for the observed variability might not be the same in all sources or that it can be a combination of more than one mechanism.For completeness we include a number of explanations that we have been able to reject or that are unlikely to be responsible for the extreme behaviour.Since not much can be said regarding the sources that have very few detections only in some of the bands, the discussion mostly considers the sources with the most complete data.

Rejected explanations
More data, especially multifrequency monitoring of the flares, are absolutely necessary to narrow down the possible explanations, however, based on the current data some scenarios can already be ruled out.These alternatives alone cannot explain the observed properties of our sources, but we cannot totally discard their presence in them.

Normal relativistic jets
Based on the results in this paper and in Berton et al. (2020b) it is obvious that the sources in our sample do not host persistent, continuously visible relativistic jets similar to those seen in other jetted NLS1s or any other class of jetted AGN.Several jetted NLS1s exhibit 37 GHz behaviour similar to the sources studied in this paper (Lähteenmäki et al. 2017), and all of them also show core or core-jet structures in mas-scale VLBI observations (e.g., Doi et   densities of the previously studied jetted NLS1s vary from a few mJy to hundreds of mJy, and thus are at a level that should have been easily detectable in our VLBA observations.However, the nondetections of these sources either imply that the radio core is very faint, < 0.5 mJy, or possibly absorbed (see Sect 5.2.4 and 5.3.2).We did not expect to be able to resolve the possible jet with the JVLA -except perhaps in the highest-frequency bands -since the flaring behaviour implies that we are seeing these sources at quite small angles.However, our initial assumption, again based on the observations of other jetted NLS1s, was that these sources would show flat or inverted spectra toward higher frequencies.Only one of our sources, J1228+5017, shows a radio spectrum that can be deemed flat, and none of the detected sources show any hints of an inverted spectrum in the JVLA observations.Regarding the nondetected sources, from these results we can only infer that their spectra do not turn inverted toward higher frequencies.
With these combined results we are able to reliably rule out the possibility that the variability in our sources is due to flares in a relativistic jet similar to those in other jetted NLS1s or AGN.This does not necessarily mean that the jet is absent, but in the low state it seems to be undetectable, implying that there must be also other contributors to the observed behaviour.

Kinematically young jets
These results also rule out one of our early hypotheses, that these sources would be kinematically young and have considerably high radio spectrum turnover frequencies due to that (O'Dea & Saikia 2021).The 37 GHz behaviour could be explained as radio flares superimposed on a convex radio spectrum of a peaked source (Tornikoski et al. 2001;Torniainen et al. 2005;Tornikoski et al. 2009).Obviously this is not the case, as we do not see any signs of spectra resembling those of peaked sources.Also the long-term temporal behaviour disagrees with this scenario, since several of these sources have been detectable at 37 GHz at the same flux density level for the past ∼ten years, ever since the observations first started.In case of a kinematically young source, the turnover frequency is expected to decrease very fast during the early stages of its life, staying above > 40 GHz only for 6-20 years (Berton et al. 2020b) -the kind of evolution we should be able to recognise at 37 GHz, and also at 15 GHz, as increasing or decreasing detectability, or as long-term permanent changes in the flux density levels.There are a few sources that have not been detected during the past few years even when they have been observed regularly (J1232+4957 and J1510+5547), which indicates temporal changes in these sources.Even in these cases kinematically young jets seem improbable since the evolution is not so fast that we would not have been able to detect a convex spectrum at lower frequencies with the JVLA.It should be noted that whereas kinematically young jets with SSA cannot explain the behaviour of our sources, it does not mean that the jets in these sources could not be young.

Fast radio bursts
The seemingly sporadic detectability, implying very short timescales, raised the question of whether this phenomenon could be related to fast radio bursts (FRB).FRBs are short, sub-second duration broadband Jy-level pulses of extragalactic origin (for a recent review, see Petroff et al. 2022).Several repeating FRBs have been found, and in principle they could fall into the MRO beam during an observation.In practise, it is very unlikely that such an event could account for the detections of these sources: first, the moderately long 1600-1800 s integration time used at MRO would average out even a Jy-level, sub-second pulse to an undetactable level, and second, FRBs have very steep spectra with an average spectral index of -1.5 (Macquart et al. 2019), making them fainter and even harder to detect at high radio frequencies.

Tidal disruption events
Tidal disruption events (TDE) occur when a star passes by too close to a supermassive black hole and gets disintegrated.In some extreme cases these events can result in launching of (mildly) relativistic jets, reaching luminosities around 10 42 erg s −1 , and therefore possibly bright enough to explain our 37 GHz detections (Alexander et al. 2020, and references therein).However, the timescales of TDEs are in the range of tens to hundreds of days and thus not compatible with the behaviour of our sources.Furthermore, so far a TDE has never been observed twice in the same source, and thus it seems extremely improbable that repeated detections over ten years could be due to TDEs.There are some records of partial TDEs (Campana et al. 2015) when the whole star does not get destroyed but continues to orbit the black hole, causing small TDEs once per orbit.Whereas partial TDEs could be responsible for repeated radio flares, they are unlikely to produce variability at a timescale of days.

Unlikely explanations
In the following we discuss some alternatives that are unlikely, but cannot be totally ruled out yet, or are not able to explain our sources on their own, but might contribute to the observed properties.

Observational effects
Interestingly, it seems that in all cases an inverted spectrum or a high state is seen only in single-dish observations, whereas interferometric observations show a barely flat or a steep spectrum, if the source is detected at all.This raises the question of whether the difference could be explained by contamination by nearby compact sources that the larger beams of the single-dish telescopes pick up, or by emission resolved-out with radio arrays.The first explanation -different beam sizes -can be ruled out since based on the JVLA images mapping the OVRO beam there are no other strong radio sources close to any of our targets, and thus even the largest beams (MRO and OVRO) should not suffer from confusion.
On the other hand, resolved-out emission can contribute to the discrepancy, but not explain all of it.In A-configuration the largest angular scales that the JVLA can see are approximately 5.3, 3.6, 2.4, 1.6, and 1.2 arcsec in X, Ku, K, Ka, and Q bands, respectively.In the worst case scenario, the lowest- source in Q band, this translates to 1.70 kpc.It is obvious that emission at these scales cannot explain the variability timescales seen in our sources.There can be a contribution from the resolved-out emission, but, for example, at 37 GHz based on the MRO detection threshold, it cannot exceed ∼200-300 mJy, otherwise we would be able to detect these sources much more frequently.Similarly, OVRO, with a beam of the same size as MRO, gives an upper limit of ∼10 mJy for the 15 GHz resolved-out emission.Since there are no emission sources at kpc-scale that can produce such an inverted spectrum between 15 and 37 GHz, it is reasonable to assume that the real 37 GHz flux density is at a similar or lower level than the 15 GHz flux density, suggesting that extreme variability is still present.
In addition, based on the preliminary results of our JVLA monitoring campaign of J1522+3934 using the B-configuration in X and K bands (VLA/23A-061, PI Berton), the beam size does not have a significant impact on the flux density.In the B-configuration the beam is about three times larger than in the A-configuration in both bands, and also the largest detectable angular scales -17 arcsec in X and 7.9 arcsec in K, 24.1 and 11.2 kpc at the redshift of J1522+3934, respectively -are significantly more extended than in A-configuration.However, the observed flux densities in A-and B-configurations are the same within the errors, further supporting that any resolved-out emission is not able to explain the difference.

Precessing jet
One alternative to explain variability in AGN is the precession in the jets (e.g.Kudryavtseva et al. 2011), leading to changes in the viewing angle and thus in the strength of relativistic boosting.Precession can be caused by a tilted accretion disk via different mechanisms, such as the radiation-driven warping instability (Pringle 1996) or the Bardeen-Petterson effect (Bardeen & Petterson 1975) due to Lense-Thirring precession (Thirring 1918).Precession can also be observed in binary supermassive black hole systems (Begelman et al. 1980).However, in all these cases the expected, and so far observed, precession period is of the order of years (e.g.Kudryavtseva et al. 2011;Liska et al. 2018;Horton et al. 2020), rather than days as in our case.It is therefore unlikely that precession on its own could explain the properties of these sources.

Intermittent activity
The lack of detectable jets in these NLS1s might indicate a kinematically young age -that was already discussed in Sect.5.1.2-or intermittent activity.Intermittent activity due to radiation pressure instabilities in the accretion disk was evoked to explain the excessive number of kinematically young radio AGN, such as GPS sources, and especially their subclass of compact symmetric objects (CSOs Czerny et al. 2009).For a black hole with a mass of 10 8  ⊙ the duration of the activity phases is estimated to be 10 3 -10 4 years, and the breaks between them 10 4 -10 6 years.For lower black hole mass sources, such as NLS1s, these timescales are shorter, but certainly not short enough to explain the variability we are observing.
Also 3D general-relativistic magnetohydrodynamic (GRMHD) simulations have yielded similar results; Lalakos et al. (2022) find that before establishing stable, powerful relativistic jets an AGN can go through several cycles of intermittent activity, with the jets turning on and off and drastically changing direction.This leads to an X-shaped radio morphology seen in 5-10 per cent of radio galaxies, and, naturally, considerable variability.Using the results in Lalakos et al. (2022) we can estimate that the launch-to-quench timescale for a black hole with a mass of 10 7  ⊙ is 10-100 years, and the jets re-emerge after 100-1000 years.The timescale is too long for our sources, but it suggests that in lower black hole mass AGN we could be able to follow, at human timescales, the chain of events from the initial launch of the jets until they are quenched by the infalling gas.As low black hole mass jetted sources NLS1s could be an optimal target for these kind of studies.
Shorter timescale intermittency can manifest itself as a result of changing injection rate of plasma into the jet base / jet (e.g.Lohfink et al. 2013;Fedorova & Del Popolo 2023).Between these events the jet can be totally absent or very weak, possibly explaining the low state of our sources.What remains unclear is whether these kind of events can account for the required short timescales and high variability amplitudes, and how these events manifest themselves in the radio regime.The classical viscous and thermal timescales associated with an accretion disk around a black hole with a mass of ∼10 7  ⊙ are too long to explain the variability, whereas the magnetic timescale dominating the inner parts of the disk can be considerably shorter (Livio et al. 2003;King et al. 2004).The magnetic timescale is the time on which the poloidal magnetic fields in different parts of the disk can spontaneously align, possibly changing the dissipation in the disk and its coupling to the jet.Local changes in the magnetic field alignment can cause small-amplitude flickering at very short timescales, whereas large-amplitude events, where the magnetic field is aligned in a considerable fraction of the disk, are more rare.Thus this kind of intermittency could possibly explain either the short timescales or the high amplitudes, but not both.
It is worth noting that even if intermittent activity would not be the culprit in this case, we do see signs of that among these sources.Assuming that we are now observing the jets in our sources at small angles as indicated by the variability, it is evident from the misalignment between the radio emission and the host galaxy in Fig. 5, panel c) in Järvelä et al. (2021) that J1522+3934 has experienced an earlier activity period.However, the projected size of the structure is almost 20 kpc, well beyond the host galaxy, implying that the activity period has been longer than what would be expected in the aforementioned scenarios.Based on the current data we also cannot determine whether the jets turned off or just changed direction.

Pure FFA
A possible way to explain the flares is to assume that the underlying radio emission of the relativistic jet is totally free-free absorbed by ionised gas in the low state, and would only occasionally break through the absorbing screen due to intrinsic flaring, or due to very fast drops in the absorption (see Sect. 5.3.2).By solving the transfer equation, it is possible to prove that such a scenario is not impossible, as it does not require an unreasonable amount of gas.Let us assume that the radiation produced by the jet is free-free absorbed as follows: where   is the optical depth,  ,0 is the radiation produced by the jet, and   is the radiation we observe after it has crossed the ionised gas.
For simplicity, let us do our calculation at 10 GHz, and assume that the jet emission is not detected.The detection threshold of the JVLA for our observations in the X band is 10 Jy, so we can assume an upper limit for the observed flux density of 30 Jy.Let us also assume that the jet has an underlying flat spectrum, and that the unabsorbed flux density at 10 GHz is 1 Jy.Using the previous equation, we can obtain an optical depth   ∼ 10.The optical depth of the ionised gas cloud depends on the absorption coefficient  ff  , following where  is the size of the absorbing cloud.The free-free absorption coefficient is where   is the electron number density,   is the number density of the ions,   the electron temperature,  the atomic number, and  ff is the Gaunt factor.Assuming hydrogen gas (  =   ), and using the approximation of the Gaunt factor between 0.3 and 30 GHz, the coefficient becomes If we integrate this assuming that the cloud has a uniform density and temperature, the optical depth becomes Inverting this equation, we can derive Since we now know that  ff  ∼ 10, we can try to calculate the size of the absorbing clouds by assuming different values of electron density and temperature, at the frequency of 10 GHz.For   = 10 4 cm −3 and   = 10 4 K, which are rather typical values, we obtain  = 38 pc.For a higher density, possibly similar to the conditions of a shock, of   = 10 5 cm −3 and   = 10 5 K, the size decreases to  = 8.6 pc.Such size is comparable to that of the Orion Nebula.Finally, if   = 10 5 cm −3 and   = 10 4 K, the resulting  = 0.38 pc, which is too small for a star forming region, but may be closer to the expectations of a region of gas ionised via shock by the jet itself.Due to the  2.1 dependence, the required size of the ionised cloud increases at higher frequencies.For example, at 50 GHz it would need to be ∼30 times larger to effectively absorb all the emission.This would imply sizes of hundreds of parsecs, unlikely ionised by the AGN, but of a characteristic size for a star forming region (e.g., Congiu et al. 2023).Lower densities and temperatures instead require unreasonably large sizes.For instance,   = 10 3 cm −3 and   = 10 4 K lead to  = 3.8 kpc, which is not realistic since this requires a uniform distribution of ionised gas as large as a small galaxy.
Even if the previous considerations show that this scenario is feasible, there are some issues that we cannot ignore.First of all, in this scenario in the low state the jet emission needs to be totally absorbed -otherwise we would see an inverted spectrum -thus the JVLA radio emission needs to originate outside the absorbed region.Were the absorption due to a star forming region, it could as well be the source of the faint low-state radio emission.As the star forming region cannot explain the variability, it would have to be intrinsic to the jet that would occasionally get bright enough to break through the FFA screen.However, assuming that the underlying relativistic jet is similar to those in other jetted NLS1s, we would assume the timescales to be comparable too, which is not the case.
Another way of producing the observed flares is by means of a variable optical depth, which in turn requires either fast moving clouds (see Sect. 5.3.2) or a rapid propagation of the jet throughout an interstellar medium with variable density and temperature.

Geometrical effects
The changes in the Doppler factor due to circumstances internal or external to the jet have been evoked to explain large-amplitude flares in AGN.Such circumstances could be result of changes in the orientation of the jet, or parts of it, or due to the jet substructure, such as a helical magnetic field (Villata & Raiteri 1999;Mignone et al. 2010;Raiteri et al. 2017Raiteri et al. , 2021)).This variability is characterised by achromatic frequency behaviour in the affected bands.
For example, an FSRQ CTA 102 has shown in the optical a somewhat similar behaviour to what we see in our sources in radio (Raiteri et al. 2017).The source increased its optical magnitude by six magnitudes, but in comparison the other frequencies were almost unaffected by the flare.In our case the flare seems to predominantly affect the radio emission and not other wavelengths (Romano et al. 2023).Raiteri et al. (2017) suggested that the variability in CTA 102 was caused by changes in the viewing angle due to peculiar jet geometry.If this is the case, we are observing different regions of the jet at different angles.In our sources only the radio emission would be seen at small viewing angle, experiencing stronger relativistic boosting due to the higher Doppler factor.This scenario could be consistent with what is seen in several simulations.Jets propagating in dense ISM cannot proceed in a straight line but tend to wiggle around the least resistance path (Wagner et al. 2012).
To estimate the feasibility of this scenario, we can estimate the level of change in the Doppler factor required to explain the extreme variability we observe in our sources.We assume the unbeamed flux density,  0 , to be at the level of the JVLA values, and the beamed,  obs , to be close to the MRO detections.The emission is boosted as: where  = 2 −  for a continuous jet stream, and  = 3 −  for a transient emission region, such as a blob or a knot in the jet.We assume the jet spectral index to be  = 0.A few different cases of the unbeamed and beamed flux densities are shown in Table 22.The Doppler factors in case of a continuous stream are very high, but more reasonable in the case of a transient emission region in the jet.We calculated the required change in the viewing angle resulting in the estimated changes in the Doppler factor (Table 22).We did the calculation with two different Lorentz factors (Γ) characteristic for jetted NLS1s: 10 and 20 (Abdo et al. 2009).In case of the continuous stream, when  = 2, Γ = 10 is not high enough to reach the Doppler factors shown in Table 22, and even Γ = 20 yields results only in case of  = 50 (Δ = 18°), thus we list the viewing angle changes only for the  = 3 case.The required changes are not unreasonable, for example, in Raiteri et al. (2017) the viewing angle change is ∼9 °.
However, in their case the timescale of the change is of the order of  several weeks, whereas in our case it is of the order of days.Also other issues remain, as discussed below.This hypothesis requires a relativistic jet to be present, but we do not see any clear signs of this in any of our sources.In the first order approximation, in this scenario either the jet needs to change direction and consequently its Doppler factor, or new components, possibly with higher Lorentz factors, would need to be ejected.Also other factors, for example, temporal variability in the physical conditions of the jet -such as the magnetic field, and the density and energy distribution of the relativistic particles -may contribute, but their impact can be expected to be less significant.
If the changes are due to the re-orientation of some parts of the jet it is hard to explain why we observe the flaring behaviour only in radio.This might require the same part of the jet to consistently change its orientation, which does not seem likely.In this case the variability should be achromatic, which is something we cannot yet study with the current data.If the flares are due to new components ejected, we would expect to see the underlying jet also when it is not flaring, since it should be relativistically boosted also between flares unless the blobs have considerably higher Doppler factors than the continuous stream.In both these scenarios the emission comes from the whole jet and therefore require the emitting region to be very close, within the innermost parsec, to the black hole, to be able to match the estimated timescales.An alternative way of producing drastic changes in the Doppler factor only in some parts of the jet is magnetic reconnection, which will be discussed in Sec.5.3.3.
Another geometrical effect in relativistic jets that causes changes in the observed flux density is due to large-scale, ordered helical magnetic fields.If the jet is magnetically dominated, the magnetic field can drive helical streams within the jet.These streams can experience differential Doppler boosting along the jet when on one side of the helix the radiation gets relativistically boosted and on the other side it gets diminished (Steffen 1997;Clausen-Brown et al. 2011;Gabuzda 2018).
In case of a continuous stream we should be able to see the jet at all times, which is not the case, so we can assume that in this scenario the flares are caused by a blob moving in the jet, thus  = 3 in Eq. 11.For simplicity and to maximise the strength of the effect, let us assume a helical magnetic field seen exactly at the helix angle.Assuming constant  the changes in the flux density only depend on the Doppler factor whose value depends on the angle between the helical stream within the jet and the line of sight as: where  = /, and  is the angle compared to the line of sight.In our scenario  has a minimum of 0°.Let us estimate the radius of the jet in case of the longest -folding timescale in a MRO-detected flare from Table 21; the 2017 decaying flare of J1509+6137.The flux density decreased from 970 mJy to 610 mJy in 5.97 days in our reference frame, thus in 5.97 days ×  in the source frame.Based on Table 11 let us assume that  0 = 0.1 mJy, and that  obs,max = 970 mJy, which happens when  = 0°.Using Eq. 11 we can estimate that the required Doppler factor at the maximum flux density is  max = 21.3, and at  = 610 mJy it is  610 mJy = 18.3.Using Eq. 12 with  max = 21.3 and  = 0°we get  = 0.996.Assuming  stays constant, and using  610 mJy = 18.3, we can solve for  610 mJy = 2.41°.The radius of the jet can be solved from where  is the radius of the jet, and  is the distance the blob has travelled along the arc of the outer edge of the jet.In this case  = 5.97 days ×  max × 0.996  = 3.28 × 10 15 m, and  = 7.80 × 10 16 m = 2.53 pc.This is the least extreme case, and in other flares where the changes were faster also the radius of the jet would have to be smaller to account for the variability.In cases when  min ≠ 0,  would have to be larger to result in the same , and  would have to be smaller than in the  min = 0 case.Based on other AGN, jet diameters of a few parsecs are measured at projected distances from ∼1 to 10 pc from the AGN core (Kovalev et al. 2020), and thus most likely outside the BLR.This brings us to the same question again: where is the jet when it is not flaring?Though it should be noted that we estimated the radius assuming the most favourable conditions for Doppler factor changes, thus it is likely that in reality the radii should be smaller, but by how much is unclear.

Jet -cloud/star interaction
Shocks in the interaction region of a jet and ISM can efficiently accelerate the electrons and thus increase the observed flux density of the jet (e.g.Fraix-Burnet 1992).In the case that the ISM consists of clumpy clouds only, parts of the jet might come in contact with them resulting in regions of enhanced emission that are smaller than the radius of the jet (Gómez et al. 2000).Particularly relevant in our case is the possible interaction between the jet base and BLR clouds or stars (Araudo et al. 2010;del Palacio et al. 2019;Bosch-Ramon et al. 2012).According to simulations the timescales of these events are of the same order as the estimated timescales of our sources, that is, from less than a day to a few days, and they can considerably increase the luminosity of the source.Whereas the timescales fit our observations, a BLR cloud or a massive star entering the jet is expected to produce a flare that should be observable over the whole electromagnetic spectrum, which is behaviour not consistently observed in our sources.On the other hand, based on Fermi data, there does not seem to be strong evidence pointing at the BLR photons interacting with the jet since most blazars do not show the expected high-energy cut-off (Costamante et al. 2018).However, this result can be explained if the main gamma-ray-emitting region in AGN is outside the BLR and swamps the gamma-ray emission originating inside the BLR.As a result, jet -cloud/star interaction can still cause flares observable in lower energies, for example, in the optical and radio regimes (Romero et al. 2000(Romero et al. , 2002)).
The issue of the missing jet also remains with this explanation.Although if the jet is small and embedded in the BLR clouds, also FFA could play a role in this scenario.Furthermore, since no dedicated simulations exist, it is unclear what the temporal evolution of these flares in radio is.More detailed simulations will be required to estimate if this hypothesis could provide a feasible explanation for the extreme variability of our sources.

Relativistic jet and FFA with moving clouds
In this scenario the starting point is similar to that in Sect.5.2.4,but the region of ionised gas is not uniform and stationary, but consists of moving ionised gas clouds.The AGN would be totally free-free absorbed most of the time and the flares take place when the nucleus is temporarily revealed.In other words, the behaviour we observe would be caused by a combination of obscuration and geometry, and not by an intrinsic change in the jet activity.Some support for this hypothesis was found in J1641+3454 in which no absorption was detected in X-rays just after a flare when the nucleus probably was exposed, but a possible warm absorber is seen in the X-ray spectrum when the source is in a low state (Lähteenmäki et al. in prep.)In this scenario the timescale only depends on the size of the gap in the clouds, its distance from the radio emitting source, and its orbital velocity, so the timescales can be arbitrarily short.
In this hypothesis the covering medium would most likely be ionised BLR clouds that are considerably denser and smaller than ISM clouds.The BLR clouds can be as dense as   = 10 11 cm −3 (Ferland et al. 1992) with sizes around 100-400 solar radii and thus easily able to absorb bright radio emission even at high frequencies.The covering factor of the BLR in optical is believed to be around 10-50 per cent, but reaching ∼100 per cent towards certain directions (Gaskell 2009).
However, open questions remain also in this case.This scenario requires that the jets of these sources are kinematically very young and still within the BLR, and also that their advancement is hindered enough by the BLR so that they have stayed within the BLR our whole observing period, about 10 years.Assuming a BLR outer radius of 0.1 pc, the jet propagation speed would have to be ≲ 0.03 c for this hypothesis to be viable.Wagner et al. (2012) report a jet propagation speed of 0.003-0.16c in the presence of clouds impeding its progress.Thus a slow jet could stay within the BLR up to a hundred years, easily enough for our case (Kino et al. 2021;Savolainen et al. 2023).

Magnetic reconnection
Magnetic reconnection in the jet or in the black hole magnetosphere has been evoked to explain fast variability in AGN, especially at GeV and TeV energies (e.g., de Gouveia Dal Pino et al. 2010;Giannios 2013;Kadowaki et al. 2015;Shukla & Mannheim 2020).It can account for high-amplitude variability at timescales from minutes to days.If magnetic reconnection were to take place in the jet in the form of so-called jets-in-jets or minĳets (e.g., Ghisellini & Tavecchio 2008;Giannios et al. 2009;Nalewajko et al. 2011), the jet would need to be heavily absorbed, since it remains undetected, and it would likely still be within the BLR.Proof of classic gamma-ray flares happening inside the BLR does exist (Vovk & Neronov 2013;Liao & Bai 2015), and also signs of gamma-ray pair attenuation have been found (Poutanen & Stern 2010), further suggesting that flares can happen inside the BLR.However, the research so far has concentrated on the high-energy characteristics of minĳets, and the production of radio emission and flares in the context of magnetic reconnection in the jet has not been studied, thus it is unclear whether this scenario could result in the behaviour we see in our NLS1s.
An alternative for the magnetic reconnection in the jet is the magnetic reconnection in the black hole magnetosphere (e.g., de Gouveia Dal Pino et al. 2010;Kadowaki et al. 2015;Ripperda et al. 2022;Kimura et al. 2022).The advantage of this explanation is that it does not require the presence of a permanent relativistic jet.The emission characteristics of these kind of events have been studied utilising general-relativistic magnetohydrodynamic simulations (Ripperda et al. 2022) and also development of the theoretical framework has been started (Kimura et al. 2022), but we still lack any direct evidence of this.de Gouveia Dal Pino et al. (2010) and Kadowaki et al. (2015) argue that the radio and gamma-ray emission in lowluminosity AGN can be explained with magnetic reconnection in the black hole magnetosphere, whereas blazars also require a significant contribution from the relativistic jet.Based on their model, an effectively accreting black hole with a mass of 10 7 M ⊙ and turbulenceinduced fast reconnection can show magnetic reconnection power spanning from 10 39 to 10 43 erg s −1 and thus likely enough to explain the flares in our sources.

Implications
It is evident that more data, especially simultaneous multifrequency observations of the flaring state, are required to determine the origin of the extreme variability seen in these NLS1s.Already based on the current data, the most common causes of radio variability in AGN can be ruled out, or they would require considerable fine-tuning.The strictest requirements come from the variability timescales, especially coupled with the extremely high, three to four orders of magnitude, amplitude of the variability.The timescales are extraordinarily short and therefore require a compact, milliparsec-scale, emitting region, or, alternatively, a peculiar interplay between the jet and the BLR clouds.Whereas intrinsic variability mechanisms allowing very short timescales and high amplitudes exist, most of them are still very little studied or only based on simulations or theoretical work, and lack observational evidence.To determine if any of them could explain the behaviour of our sources, more detailed theoretical framework, possibly dedicated simulations, and especially targeted observations will be needed.It should be kept in mind that we cannot exclude the possibility that we are seeing a new type of variability either.Either way, catching flares in these sources will be challenging due to the short timescales and sporadic activity, but considering that these NLS1s exhibit one of the most extreme radio variability seen in AGN so far, they do deserve our full attention.
One of the most interesting aspects of the discovery of these sources is that they were found among two very differently selected samples whose final detection percentage at MRO turned out to be very high at 12 per cent -eight sources out of 66 were detected.Whether our selection criteria actually helped us select NLS1s exhibiting this behaviour or if it was pure coincidence remains unclear.Observations of other NLS1 samples selected using similar and different selection criteria will be needed to estimate the impact of the selection effects.
Either way, detecting >10 per cent of a presumably mostly radiosilent NLS1 sample is extraordinary and raises the question whether this variability phenomenon is characteristic of NLS1s or possibly early-stage AGN, or if similar sources are hiding also among radioweak AGN of other classes.For obvious reasons radio-weak AGN have not been a target of extensive high radio frequency monitoring campaigns and we therefore know very little about their behaviour in that regime.It is possible that also strong radio AGN exhibit this kind of behaviour, but that it is swamped by other sources of radio emission and thus has remained undetected.Investigating in which kind of sources this phenomenon can be seen can help us to determine the cause of the variability.Being able to identify any common properties these sources have will also help us to find more of them.
Whether this kind of variability is limited to early-stage AGN or if it is a more common phenomenon has implications on our current understanding of AGN.These sources clearly represent an unknown population of AGN, that has gone unnoticed so far.If they are a new type of jetted AGN or something else entirely, is unclear, as is their evolution and relation to other classes of AGN.Furthermore, we do not know how common they are or if they are characteristic to the local Universe, or also exist at higher redshifts.Further studies are also required to estimate which kind of a role they play in AGN feedback over the cosmic time.

CONCLUSIONS
In this paper we investigated the origin of the extreme radio variability seen in seven NLS1s using the JVLA, the VLBA, MRO, and OVRO observations.These extraordinary sources defy an easy explanation, but the new data presented in this paper allowed us to rule out some alternatives and set additional constraints on the possible explanations.Our main conclusions are: • The behaviour of these sources is hard to explain with the usual variability mechanisms in AGN -instead a more complex scenario or possibly a new type of physical mechanism to produce variability is required.
• The amplitude of the variability -three to four orders of magnitude -seen in these sources is unprecedented, but it remains unclear whether it is intrinsic to the source, or caused by some external circumstances.
• The variability timescales indicate that if the variability is intrinsic the emitting region needs to be milliparsec in size.This implies that the emission originates close to the black hole, clearly inside the BLR, or from limited, confined regions in the jet.
• The high detection percentage among the original sample, that were not expected to be strong radio emitters, implies that this kind of sources could be quite common, but so far our understanding of this new population of AGN is very limited.
Revealing the nature of these peculiar sources is of utmost importance as they might be the first representatives of a new type of AGN variability, and/or a new class of jetted AGN entirely.In the future, an increase in the sample size will be essential to explore this new population.Their short timescales and sporadic activity pose an observational challenge, also given how diverse their behaviour is at different frequencies.High-cadence multifrequency radio monitoring with an instrument sensitive enough to detect also the rising and decaying parts of the flares will be essential to better characterise their variability, and set additional constraints to the different hypotheses concerning these sources.Furthermore, given the small spatial scales implied by the variability timescales, many of the upcoming telescopes and instruments currently under development, such as the next generation VLA (ngVLA) in radio, the Multi-Conjugated Adaptive Optics (MCAO) Assisted Visible Imager and Spectrograph (MAVIS) and the High Angular Resolution Monolithic Optical and Near-infrared Integral field spectrograph (HARMONI) in the optical/near-infrared, and Athena in X-rays, will be crucial to study the spatial properties and evolution of these remarkable sources.

A3 MRO
Measurements that are considered to be of poor quality, for example, due to unfavourable weather conditions, are discarded semiautomatically.Additionally faint detections are checked manually in the final data reduction stage.Bad weather conditions or other environmental effects are taken into account and also, for example, conspicuous but rare flux density increases caused by aircraft in the telescope beam.The general flux density levels are checked to be consistent with adjacent measurements (i.e.other sources observed before and after the target source).In addition, we checked if the observations of these sources could be contaminated by a bright radio source falling into the reference beam of the MRO telescope.Using LOFAR, FIRST, and VLASS data we concluded that whereas there are a few moderately bright sources with flux densities around a few hundred mJy at low radio frequencies that could be in the reference beam, all of them have steep spectra, and it is thus very unlikely that any of them could affect these observations.Due to the fairly high detection limit of the telescope (i.e.approx.200 mJy in optimal conditions, which is more than adequate for the bright AGN monitoring programmes conducted at MRO), we typically only see the highest tips of the flares in faint sources, whereas most of the lower level activity remains below the detection threshold (e.g., Acciari et al. 2014).This is also seen in the upper limits, the level of which can drastically change in even a short time due to compromised weather conditions that can also significantly raise the detection limit.The undetected source could be actually fainter due to variability or the observing conditions could be worse (or both), and it is therefore not detected.The upper limits describe the largest flux density the source could have in the current conditions, but still remain below the 4 detection limit and cannot therefore be used for data analysis.Occasionally there are several high upper limit values clustered around detections, which indicate that the source could indeed be active, but for some reason it does not exceed the current detection limit, for example, due to weather.However, it has been shown that the high activity periods of NLS1 sources detected at MRO correspond to flare peaks in OVRO data (Lähteenmäki et al. 2017), confirming that at least most of the time the two telescopes are catching the same events.

A4 OVRO
The AGN monitoring sample at OVRO mostly consists of bright blazar-type objects, with the majority having a mean flux density > 60 mJy (Richards et al. 2014).Therefore in the schedules, each observation consists of four on-on integrations, each 8 s long, resulting in a total integration time of 32 s.Given the small number of on-on integrations, it is possible that atmospheric fluctuations or pointing errors result in outliers in the light curves (Richards et al. 2014).Moreover, the number of observations performed each day is large (up to 500) so that it is possible that some, apparently high signal-to-noise ratio observations, occur purely due to random fluctuations.The data are processed with an automated pipeline, where manual editing is done to flag obviously bad periods of data, and data are automatically flagged based on large changes within the four on-on integrations and other diagnostics (see Richards et al. 2011, for details).However, individual data points are not typically manually inspected, as for the variability analysis of bright blazars, the effect of outliers in the data is small (Richards et al. 2014).
Because of the faintness of the NLS1 targets, we have done additional manual checks to inspect the quality of the detections, which we describe here.We note that in all the cases described below, the flux density of the spurious detection has been less than 20 mJy and mostly < 10 mJy, or the uncertainty has been larger than usual so that similar observations in our blazar light curves would not be as problematic.
The receiver records both right-and left-hand circular polarisation separately with the final observation being a weighted average of the two.We can thus inspect the two values separately to verify that the source has been detected in both polarisations (here we assume that the circular polarisation of the objects is negligible, as is the case for most blazars at 15 GHz (e.g., Homan & Lister 2006).This made us reject two spurious detections in J1232+4957.Additionally, we have inspected observations of other nearby sources to see whether there are data that have been automatically flagged in the pipeline close to the observation of the NLS1, indicative of potentially poor observing conditions.This resulted in the rejection of one spurious detection in J1641+3454.
In October 2021 we also changed the observing strategy of these NLS1 targets so that they are observed twice in a row in the sched- ules.This way we can see whether short-term atmospheric effects or bad conditions have resulted in spurious detections if the two consecutive observations show a large difference, as we would not expect large changes on a time scale of ∼ 1 min.This resulted in the rejection of single spurious detections in J1029+5556, J1509+6137, J1510+5547, and J1522+3934, all of which had moderate S/N values of ∼4-9.
The remaining detections in the paper either show detections in two consecutive observations (J1522+3934) or consistent flux densities in the right-and left-hand circular polarisation and no apparent problems with nearby targets (J1029+5556, J1522+3934, J1641+3454).However, we cannot fully exclude additional, unknown effects in the observations before October 2021 when the sources were observed only once in a row, especially in J1029+5556 and J1641+3454 that do not show any other detections in the OVRO light curves.J1522+3934 on the other hand seems more reliable because of its multiple detections.

APPENDIX B: RADIO MAPS
The JVLA radio maps with likely detections are shown here.This includes X, Ku, K, Ka, and Q band maps of J1228+5017 (Figs.

Figure 1 .
Figure 1.Non-simultaneous radio spectrum of J1029+5556.Symbols and colours explained in the figure.Filled symbols denote integrated flux densities and empty symbols mark peak flux densities, except empty red symbols with downward arrows that are used for upper limits.VLA 1 data from Berton et al. (2020b) and VLA 2 data from this paper.

Figure 2 .
Figure 2. Non-simultaneous radio spectrum of J1228+5017.Symbols and colours explained in the figure.Filled symbols denote integrated flux densities and empty symbols mark peak flux densities, except empty red symbols with downward arrows that are used for upper limits.VLA 1 data from Berton et al. (2020b) and VLA 2 data from this paper.

Figure 3 .
Figure 3. Non-simultaneous radio spectrum of J1232+4957.Symbols and colours explained in the figure.Filled symbols denote integrated flux densities and empty symbols mark peak flux densities, except empty red symbols with downward arrows that are used for upper limits.VLA 1 data from Berton et al. (2020b) and VLA 2 data from this paper.

Figure 4 .Figure 5 .
Figure 4. Non-simultaneous radio spectrum of J1509+6137.Symbols and colours explained in the figure.Filled symbols denote integrated flux densities and empty symbols mark peak flux densities, except empty red symbols with downward arrows that are used for upper limits.VLA 1 data from Berton et al. (2020b) and VLA 2 data from this paper.

Figure 6 .
Figure 6.Non-simultaneous radio spectrum of J1522+3934.Symbols and colours explained in the figure.Filled symbols denote integrated flux densities and empty symbols mark peak flux densities, except empty red symbols with downward arrows that are used for upper limits.VLA 1 data from Berton et al. (2020b) and VLA 2 data from this paper.

Figure 7 .
Figure 7. Non-simultaneous radio spectrum of J1641+3454.Symbols and colours explained in the figure.Filled symbols denote integrated flux densities and empty symbols mark peak flux densities, except empty red symbols with downward arrows that are used for upper limits.VLA 1 data from Berton et al. (2020b) and VLA 2 data from this paper.
B1-B5), X, and Ku band maps of J1232+4957 (Figs.B6-B7), X, Ku, K, and Ka band maps of J1522+3934 (Figs.B8-B11), and X, Ku, K, Ka, and Q band maps of J1641+3454 (Figs.B12-B16).APPENDIX C: LIGHT CURVES The light curves of our sources from the beginning of 2014 to mid-2022 are shown here.Figs.C1-C7 show light curves including lowresolution (MRO and OVRO) and high-resolution (JVLA, VLBA, and VLASS) data.Due to the strongly varying flux densities these plots are in logarithmic scale.The light curves in Figs.C8-C14 show only the MRO and OVRO data in linear scale, and include also the 4 upper limits for both observatories.

Table 2 .
Summary of the single-dish observations published here.

Table 22 .
Required Doppler factors and changes in the viewing angle of the jet.Columns: (1) unbeamed flux density; (2) beamed flux density; (3) required Doppler factor assuming a continuous jet stream; (4) required Doppler factor assuming a moving component in the jet; (5) required change in the viewing angle assuming  = 3 and Γ = 10; (6) required change in the viewing angle assuming  = 3 and Γ = 20.