Unstable Phenomena in Stable Magnetospheres: Searching for Radio Flares from Magnetic OBA Stars Using VCSS

Although the majority of hot magnetic stars have extremely stable, ∼kilogauss strength surface magnetic fields with simple topologies, a subset undergoes small-scale explosions due to centrifugal breakout. The resulting small-scale flares are typically below the sensitivity of current magnetospheric diagnostics and do not generate detectable transient signatures. However, a recently reported radio flare from the hot magnetic star CU Vir suggests that some of the most energetic events do reach detectable levels. Motivated by this, we searched for transient radio sources in the first two epochs of the VLITE Commensal Sky Survey at the positions of 761 hot magnetic stars. We report three detections. A false-association analysis shows a less-than-1% probability that the sources are imaging artifacts. We then examine the stellar parameters of the three stars to understand if they are likely to produce flares. We conclude that while at this stage, we cannot make a definitive association of the detections with the stars, the current data are consistent with the hypothesis that the flares originate in the stellar magnetospheres.


INTRODUCTION
Magnetic fields have been observed in planets, brown dwarfs, and stars of all spectral types, from ultracool dwarfs to the very hot massive stars (spectral types O, B and A, or early-type stars) (e.g.Grunhut et al. 2012;Shultz et al. 2018;Hallinan et al. 2006;Kao et al. 2018;Zarka 2004, etc.).The characteristics of stellar magnetic fields vary with spectral type, which is often attributed to the origin of the field in a given class of star.For a Sun-like star, the field is produced and sustained by a convective dynamo operating close to the surface of the star (e.g.Tobias 2002).These fields are highly dy-namic and have complex topologies.They give rise to explosive phenomena, such as flares and coronal mass ejections (CMEs).These phenomena are collectively referred to as stellar activity.Note that all stellar activity is inherently associated with some change in the magnetic configuration (reconnection of field lines).
By contrast, the magnetic OBA stars, which constitute ≳ 7% of the early-type star population1 (Grunhut et al. 2012(Grunhut et al. , 2017;;Sikora et al. 2019;Hubrig et al. 2023), are thought to have quieter surroundings.This is due to the fact that their magnetic fields are extremely stable throughout their lifetimes, and typically have simple topologies (e.g.Kochukhov et al. 2019).In most cases, the magnetic fields in OBA stars can be described as a dipole tilted to the rotation axis (Stibbs 1950;Shultz et al. 2018;Sikora et al. 2019, etc.).The interaction of these stable magnetic fields with the stellar winds results in various magnetospheric phenomena that exhibit only periodic variability due to rotational modulation arising from the misalignment between rotational and magnetic axes (e.g.Oksala et al. 2015;Leto et al. 2017).This predictability represents a fundamental difference from the random stellar activity of late-type stars.
In recent years, ideas regarding the magnetospheric dynamics in rapidly rotating and strongly magnetic hot stars have undergone revolutionary changes.Such stars have magnetospheres in which the Kepler radius, R K (the distance at which the centrifugal force due to corotation balances gravity), is smaller than the Alfvén radius, R A (the radius of the largest closed magnetic field line).The region inside R K is called the Dynamical Magnetosphere, and that between R K and R A is called the Centrifugal Magnetosphere (DM and CM respectively, Petit et al. 2013).
In the DM, the radiatively driven stellar wind material follows the magnetic field lines, and is then drawn back to the star (again following the field lines) by gravity.Inside the CM however, the outward centrifugal force is stronger than gravity, and it is the magnetic tension that keeps material from escaping the stellar magnetosphere.By considering a rigidly rotating magnetosphere, Townsend & Owocki (2005) showed that within the CM, there are certain locations along the field lines where plasma can stably accumulate.The scenario of stable plasma clouds was successfully used to explain the observed periodic variability of different magnetospheric phenomena (e.g.Oksala et al. 2015).The only limitation of this scenario was that it did not explain what balances the continuous accumulation of stellar wind plasma inside the magnetosphere, since the magnetic field has only a finite capacity to hold material against the centrifugal force.
The first explanation proposed was that once that stage is reached, a large-scale explosion takes place in which the magnetosphere opens up and all the confined plasma escapes (Townsend & Owocki 2005).This phenomenon was named centrifugal breakout (CBO) and predicted a highly variable magnetosphere surrounding these stars in which large-scale explosions occasionally occur.But no variability in the optical light curve other than that incurred due to rotational modulation was detected; nor were systematic changes in the depths of the light minima due to reductions in the magnetospheric column density (Townsend et al. 2013).
An alternate explanation invokes more complex physical processes such as diffusion of materials through the magnetosphere (Owocki & Cranmer 2018).This scenario predicts that the properties of different magnetospheric phenomena (such as the emission strength of the Hα emission) are dependent on the stellar mass-loss rate or its proxy.In 2020, by conducting comprehensive examination of Hα emission properties for stars with CM, Shultz et al. showed that this was not true.
The indifference of the magnetospheric phenomena to the mass-loss rate is, however, easily explained in the CBO scenario.This realization led to the new proposition that massive star magnetospheres are always maintained at the highest possible plasma density such that CBOs occur continuously, not as one single large-scale explosion that empties the magnetosphere but rather as a collection of small-scale eruptions occurring along different magnetic azimuths (Shultz et al. 2020;Owocki et al. 2020).Note that this study did not include late Btype or A-type stars as such stars are not usually found to be bright in Hα, presumably due to their weak stellar wind, or metallic wind (Shultz et al. 2020, and references therein).This raises the possibility that the cooler earlytype stars might differ from their hotter counterparts in terms of magnetospheric operation and phenomena, and CBO may not be relevant for late B/A-type stars.
This was, however, found to be unlikely when Leto et al. (2021) and Shultz et al. (2022) showed that the incoherent radio emission from early-type stars, including down to A0, follows the same scaling law with stellar parameters.This observation can be satisfactorily explained by invoking CBO driven magnetic reconnections (Owocki et al. 2022).
Thus, the most recent idea is that centrifugal magnetospheres of OBA stars are not quiet, but that they actually experience numerous tiny explosions leading to the ongoing escape of the magnetospheric plasma.Because we cannot resolve these objects, the signature of these continuous explosions, distributed randomly around the stellar magnetosphere, get washed out, giving the deceptive appearance of a stable magnetosphere.
Observationally, evidence for these small-scale explosions has been limited.Due to its high directivity and consequent sensitivity to changes on small spatial scales, the detection of electron cyclotron maser emission (ECME) is one possibility (Das & Chandra 2021).Alternately, it is tempting to extrapolate from solar flares, for which the occurrence frequency of the flares follows a power law with flare energies such that weaker flares are more frequent than stronger flares (see Figure 10 of Aschwanden et al. 2000).If a similar case holds for CBOs in hot star magnetospheres, we would expect that there will be some instances when a strong explosion will occur that will not get washed out post spatial averaging, and hence will be detectable as a flare.
In the past, there have been reports of observation of flares from massive stars in various wave-bands such as optical, X-ray and radio ( §2).Some of these stars were subsequently discovered to have low-mass companions which might have hosted the flare (e.g.Bouy et al. 2009;Pedersen et al. 2017).Even for the cases when no companion was detected, the flares were considered the manifestation of an otherwise invisible companion (e.g.Nazé et al. 2014).The key reason behind this was that there was no theoretical explanation for the production of flares in hot star magnetospheres.
With the new understanding provided by the CBO mechanism it is now important to revisit the idea of whether hot magnetic stars flare or not.Since one cannot rule out the possibility of a companion star or a foreground object for the individual cases, the only way to understand whether or not hot magnetic stars flare is via statistical analysis.This paper describes the first attempt to achieve that goal.
In this work, we present the first results from a search for radio flares from hot magnetic stars using the Very Large Array (VLA) Low-band Ionosphere and Transient Experiment (VLITE) Commensal Sky Survey (VCSS)2 .This paper is organized as follows.In Section 2 we review previous claims of flare detections from hot (OBA) stars.Section 3 describes our strategy for probing transient phenomena from hot magnetic stars.A description and analysis of VCSS data is given in Section 4. Section 5 presents a discussion of our results with a summary given in Section 6.

PREVIOUS CLAIMS OF DETECTION OF FLARES FROM HOT MAGNETIC STARS
X-ray flares have been reported from both magnetic (field strength ≳ 100 G) and non-magnetic/weaklymagnetic hot stars.One of the earliest observations of flares from the direction of a magnetic early-type star was made by the ROSAT survey in February 1995 for the magnetic B star σ Ori E (Groote & Schmitt 2004).From the same star, another X-ray flare was reported by Pallavicini et al. (2002) and Sanz-Forcada et al. (2004) using the XM M − N ewton telescope in March 2002.Subsequently however, a low-mass companion was discovered by Bouy et al. (2009), and the previously observed flares were attributed to the companion as the latter belongs to an established class of flare emitters.This was disputed by Mullan (2009), who considered the low possibility of occurrence of two large flares on a late-type star within a span of just seven years.Petit et al. (2012) observed the star at higher spatial resolution with the Chandra X-ray Observatory, and found the low-mass companion to have significant contribution (40%) to the overall X-ray flux observed from the system.This discovery proved that the low-mass companion produces significant amounts of X-rays and, as a result, made the claim of σ Ori E as the origin of the flare more doubtful.
In addition to σ Ori E, there are two other hot stars, both with ∼kG strength surface magnetic fields (Shultz et al. in prep.) and not known to have low-mass companions, from which X-ray flares were reported.These two stars are IQ Aur (Robrade & Schmitt 2011) and HD 47777 (Nazé et al. 2014).In both cases the presence of an undiscovered low-mass companion was speculated to be the source of the observed flares.
Apart from X-rays, the only other wave-band where a flare has been observed from the direction of a magnetic massive star is radio.Trigilio et al. (2000) first discovered that the magnetic early-type star CU Vir produces periodic radio pulses by ECME.In the same observations they also observed secondary enhancements at 1.4 GHz that could not be confirmed as persistent features.In 2021, Das & Chandra reported ultra-wideband radio observations of the same star that cover the frequency range of 0.4-4.0GHz.In addition to the expected periodic radio pulses, they also observed strong (comparable to the persistent radio pulses), highly circularly polarized enhancements at rotational phases much offset from those of the persistent radio pulses at sub-GHz frequencies.Their subsequent observations conducted after a year confirmed these enhancements to be transients.This is the first (and only) confirmed observation of radio flares from the direction of a magnetic massive star.

STRATEGY
Our aim is to understand whether or not hot magnetic stars can exhibit transient phenomena (flares).There are two ways to attempt to achieve that goal: the first is to monitor a selected sample of hot magnetic stars (targeted observations), and the other is to look for such flares in existing databases from all-sky surveys.For our science, we deem the latter to be more efficient since it will include a much larger number of hot magnetic stars than what is possible through targeted observations.
There are two wavebands in which flares from hot magnetic stars have been claimed in the past: X-ray and radio.As a first step we chose to work with radio surveys.A number of large area radio surveys are either ongoing or recently completed over a wide range of frequencies, including the VLA Sky Survey (VLASS, 3 GHz; Lacy et al. 2016), ThunderKAT survey (1.8 GHz, Fender et al. 2016) While we do not have information to prefer one frequency range over another, we consider two important aspects that have been recently reported regarding nonthermal radio emission from hot magnetic stars.The first is that Leto et al. (2021) discovered incoherent radio emission from hot magnetic stars have spectra that exhibit turn-over at ∼ 1 GHz and the flux density decreases as one goes to frequencies below 1 GHz.This is consistent with the low detection rate of Ap-Bp stars in the LoTSS survey (Hajduk et al. 2022).On the other hand, for coherent radio emission from hot magnetic stars, the sub-GHz frequency range above 200 MHz appears to be more preferred than the GHz regime (Das et al. 2020;Das & Chandra 2021;Das et al. 2022a,b).In addition, the only confirmed observation of a radio flare from the direction of a hot magnetic star was obtained at 0.7 GHz (Das & Chandra 2021) and a tentative flare at 0.4 GHz was also reported by Das & Chandra (2021).Considering these aspects, we chose to work with a survey that operates above 200 MHz and below 1 GHz.The VCSS survey meets this criterion.We give details of VCSS in the next section.
We use a catalog of magnetic hot stars assembled from an exhaustive literature compilation (Shultz et al. in prep.).The catalog includes 761 stars of spectral types O, B, and A with at least one confirmed magnetic detection.Stars on the main sequence, pre-main sequence, and post-main sequence are included.The catalog also provides effective temperatures, luminosities and rotation periods, where available.This information is used to assess whether stars with observed transient radio emission have centrifugal magnetospheres necessary for the bursts to originate from CBO-driven flares.
Note that in the event of a flare detection, our current approach is not sufficient to discern whether that flare originated at the stellar magnetosphere or a magnectically active companion, but it does provide candidates that will be objects of interest for follow-up targeted observation campaigns.

THE VCSS DATA
VLITE is a commensal instrument on the Karl G. Jansky VLA developed by the Naval Research Laboratory (NRL; Clarke et al. 2018).VLITE collects data at frequencies between 320-384 MHz from the low frequency prime focus receivers during nearly all regular VLA operations at GHz frequencies.VLITE data are recorded on up to 18 antennas, correlated, and transferred to NRL for automated pipeline processing and archiving.Removal of persistent radio frequency interference in the upper portion of the VLITE band gives an effective bandwidth ∼ 40 MHz centered on 340 MHz.
A special correlator mode was enabled to allow VLITE processing of the on-the-fly data recorded during observations for VLASS.VLASS is a multi-epoch survey observing the δ > −40 • sky three times starting from September 2017.We use VLITE data from the first two observing epochs in this work.
VLASS tessellates the sky with observing tiles covering 10 • × 4 • in right ascension (RA) and declination, respectively.Tiles are observed by scanning in RA along lines of constant declination spaced by 7.2 ′ , the half power radius of the VLA primary beam at 3 GHz.VLITE data are allowed to accumulate using a standard 2 s sampling as the antennas move through an angular distance of 1.5 • .This corresponds to ∼ 28 s of time over most of the sky for the standard VLASS slew rates.These data are then correlated and imaged using the midpoint of the motion as the phase center.The result is short "snapshots" with a smeared or elongated primary beam response.The Perley & Butler (2017) fluxes for 3C138 and 3C286 are used to set the flux scale.The calibrated and imaged VLITE data recorded during VLASS define VCSS.
Approximately 160, 000 snapshot images covering a ∼ 32, 000 deg 2 sky footprint are produced during each VCSS epoch.Snapshots are 3.5 • across with a typical sensitivity 7 − 10 mJy bm −1 that increases away from image center due to primary beam attenuation.The attenuation is about a factor of three at the edge of the field of view (FoV).
VCSS snapshots are highly overlapping due to VLITE's wide FoV (3.5 • ) relative to the declination spacing (7.2 ′ ) of the VLASS observing strips.As a result, each point on the sky is typically sampled by 40−50 snapshots spread over 80 − 100 min, although irregular sampling over several months to a year occurs along tile boundaries.
Snapshots in the first survey epoch are imaged with elliptical synthesized beams with axes ranging 7 − 35 ′′ .In the second epoch, snapshots are imaged with 20 ′′ round beams to facilitate combination into more sen-sitive mosaic images.While these are useful to survey for long timescale transients, combining snapshots imaged over hours to months suppresses short timescale emission typical of stellar flares.For this reason, we do not use the mosaic images in this work and focus our search on the VCSS snapshots.
VCSS snapshots are processed with the VLITE Database Pipeline (VDP, Polisensky et al. 2019).The Python Blob Detector and Source Finder (PyBDSF, Mohan & Rafferty 2015) is used to catalog sources with signal-to-noise ratios (SNR) greater than five, where the noise is calculated locally around the source position.
Querying the VDP database for snapshot sources within 10 ′′ (the approximate minimum VCSS resolution) of hot magnetic stars from the catalog revealed three matches.Each VCSS match appears in only one snapshot and none correspond to known radio sources.However, the VCSS survey is known to contain a large number of imaging artifacts.More than one million sources in each epoch are detected only once and unassociated with known radio sources.We therefore performed an analysis to determine the likelihood our star detections are due to false associations.

VCSS Artifacts
In Appendix A, we provide an in-depth analysis of artifacts present in VCSS snapshots.Our findings reveal the existence of two distinct artifact populations: one concentrated within an arcminute of bright sources, and another distributed uniformly in all directions.Given that the nearest neighbors for the three star detections are situated more than 500 ′′ away, we concentrate our analysis on the isotropic artifact population.
We also observe a correlation between the density of artifacts and the SNR; fewer artifacts are observed at higher SNR.Additionally, the substantial variation in the total artifact count per snapshot indicates the number of false associations within any stellar catalog will depend on the detection threshold, the distribution of stars and the specific set of snapshots they encompass.

False Detection Probability
The probability, P , of observing N events when λ are expected is given by Poisson statistics: The expected number of false associations in a single snapshot is determined by the product of the solid angle searched and the density of artifacts in the snapshot.All three VCSS detected hot magnetic stars have SNR > 5.5.Therefore, we calculate the product of the 10 ′′ search radius solid angle, Ω * , with the isotropic artifact density above a SNR of 5.5 for each snapshot, n/Ω F oV , where Ω F oV is the snapshot solid angle (9.6 deg 2 ).We sum over all stars and all snapshots they are situated in to get the expected number of false associations for the catalog: We calculate λ = 0.44 for our star catalog.The probability of three false associations when 0.44 are expected is 0.9%.We confirm this result with simulations of randomly generated star catalogs in Appendix B.

RESULTS AND DISCUSSIONS
The three stars detected by VCSS are HD 36644, HD 40759 and HD 175362.Figure 1 presents the corresponding light curves for these stars.Brightness errors include 1σ fitting uncertainties with 3% primary beam correction uncertainty added in quadrature but do not include the 15% VLITE flux scale calibration uncertainty.
Upper limits from non-detections are defined as three times the local noise.We use the local noise maps calculated by PyBDSF across the FoV corrected for primary beam attenuation.Implementing rigorous quality controls, VDP prioritizes the flagging of poor quality images, even if it results in a small fraction of quality images being flagged.In cases where our examination indicated reliable flux measurements from flagged images, their upper limits have been included in Figure 1.A discussion of each detection follows.
HD 36644 was detected on 2017 October 2 at 09:19:25 UT during a 28 s duration snapshot.The brightness was 107 ± 8 mJy bm −1 with a SNR of 5.6.Its nearest neighboring source was 918 ′′ away.Non-detections effectively constrain the emission timescale to under seven minutes.
HD 36644 lay beyond the northern border of the tile under VLASS observation on the day of the detection.The stellar position, however, falls within the coverage of nine VCSS snapshots spanning six VLASS declination strips.These strips were observed sequentially from north to south, resulting in a general trend of HD 36644 appearing farther out in the VCSS primary beam as time progressed.The larger beam attenuation at larger distances creates the upward drift in the non-detection limits in Figure 1.
For illustrative purposes, Figure 2 presents cutout images centered on the position of HD 36644 at the times indicated in the light curve.Notably, intensity noise at the stellar location increases as the distance from the snapshot center grows, primarily due to the effects of primary beam attenuation.
HD 40759 was detected with a significance of 6.8σ on 2020 October 9 at 15:21:23 UT.Similar to the previous case, this detection occurred in a 28 s snapshot, with the source positioned 1248 ′′ from its nearest neighbor.Its observed brightness was 113 ± 9 mJy bm −1 .The upper panel of Figure 3 presents the image cutout captured at the time of detection.
During the time of its detection, the star was situated off the northwest corner of the tile under observation by VLASS.Over time, observations along the strips proceeded south to north, gradually approaching the stellar location.The decreasing primary beam attenuation cre-ates the downward trend in non-detection limits within the HD 40759 light curve.While the duration of the emission remains unconstrained, the non-detections preceding the burst effectively impose an upper limit on the rise time of ∼ 5 minutes.
Unlike the other two stars, HD 175362 was detected within the tile being observed by VLASS.The temporal pattern of non-detections with values below ∼ 50 mJy bm −1 exhibits a distinctive sequence of decreasing, flattening, and then increasing.This pattern corresponds to the primary beam attenuation sampled across the north-south direction of the VLASS declination strips by snapshots aligned with the stellar RA.Additionally, the larger non-detection values occurring at all times reflect the larger primary beam attenuation of snapshots centered at large distances in the east-west direction along the observing strips.
The detection of HD 175362 occurred on 2022 February 8 at 16:43:43 UT during a 42 s snapshot.The fitted brightness was 69 ± 10 mJy bm −1 , 5.5 times the local noise level.Positioned 529 ′′ away from its nearest neighbor, this detection is flanked by non-detections before and after it, effectively constraining the emission timescale to less than ∼ 1 minute.To further illustrate this event, the lower panel of Figure 3 presents the image cutout at the moment of detection.
Although none of the stars are detected by VLASS, all VCSS detections occurred beyond the narrow VLASS observing strip coincident with each snapshot.Specifically, HD 36644 remained unobserved by VLASS until 10 days after its VCSS detection.Similarly, HD 40759 had been observed by VLASS 82 days prior to its VCSS detection.The VCSS detection of HD 175362 preceded the VLASS observation of the star's position by ∼ 20 minutes.However, the VCSS upper limit on the burst timescale is considerably shorter than the VLASS observing delay.
In the following subsections, we discuss the plausibility of these stars serving as potential sources of the observed radio transients.

Do the stars have centrifugal magnetospheres?
As described in the introduction, CBOs are the only known phenomena that could give rise to an observable transient event from a hot magnetic star without any influence of a companion.But for this scenario to be valid, the stars must have centrifugal magnetospheres (CMs).In this subsection, we examine whether the three stars harbor CMs.
The two parameters needed to examine this property are the Kepler radius R K and the Alfvén radius R A .R K can be calculated using the following equation; where, G is the universal gravitational constant, M * is the stellar mass and Ω = 2π/P rot is the stellar rotational frequency.R A , on the other hand, is a more complicated function of several stellar parameters such as the stellar magnetic field strength, rotation period, mass-loss rate, wind speed and also the misalignment between rotation and magnetic axes (obliquity β).R A increases with increasing magnetic field, rotation period, and decreases with increasing mass-loss rate.Barring the case of the very rapid rotators (P rot ≲ a day), R A decreases with increasing stellar wind-speed, and is nearly independent of β.For the very rapid rotators, R A can increase with increasing stellar wind-speed and increasing β (e.g.see §2.2 of Trigilio et al. 2004).As can be seen below, none of the three stars under consideration are very rapid rotators, so that we can ignore the dependence of R A on β and the rotation period.In that case, the equation for (4) where R * is the stellar radius, and η * is the 'magnetic confinement parameter' (Ud-Doula et al. 2008) that is a function of the surface equatorial magnetic field B eq , R * , mass-loss rate Ṁ and wind terminal speed v ∞ .Among the three stars, HD 175362 is the most wellstudied.It is also the hottest star among the three (T eff = 17.6 ± 0.4 kK, Shultz et al. 2019a).The star has a magnetic field with polar strength of ≈ 17 kG (Shultz et al. 2019b).Its rotation period was estimated as 3.67381(1) days (Bohlender et al. 1987;Shultz et al. 2018), where the uncertainty in the least significant digit is in brackets.The stellar radius and mass were estimated to be R * = 2.7 ± 0.2 R ⊙ and 5.3 ± 0.2 M ⊙ (Shultz et al. 2019b), respectively, which gives a Kepler radius of 6.5 ± 0.4 R * .This is much smaller than the reported Alfvén radius of the star which is ≈ 46 R * (Shultz et al. 2019b).Thus the star harbors a centrifugal magnetosphere, and hence may experience CBO phenomena.This is also supported by the fact that the star's non-thermal gyrosynchrotron radio luminosity follows the scaling relation with stellar parameters expected for CBO driven radio emission (Shultz et al. 2022;Owocki et al. 2022).
The second best-studied star is HD 40759, a de-facto spectroscopic triple-star system with all three components being A1-B9 stars.The effective temperature T eff of the components lies in the range 9-12 kK (Semenko et al. in prep).The magnetic component is a chemically peculiar star with a rotation period of 3.37484 days (Semenko et al. 2022).The maximum observed longitudinal magnetic field is ≈ 2 kG, with a root mean square value of 1.12 kG (Romanyuk et al. 2021;Semenko et al. 2022).Accurate mass and radius estimates for the star are not available at the moment, yet we can evaluate them roughly from the evolutionary status of the star.HD 40759 is a member of the stellar association Orion OB1, subgroup 1c, with an average age of 4.6 Myr (Semenko et al. 2022).Interpolating stellar evolutionary tracks provided by the MIST project (Dotter 2016;Choi et al. 2016), gives the mass and radius of the star as 2.5 M ⊙ and 2.1 R ⊙ .Thus, we obtain R K = 6.1R * .
From the lower limit of the surface magnetic field (2 kG), we can also estimate the Alfvén radius by assuming typical values of mass-loss rate ( Ṁ ) and wind terminal speed (v ∞ ) observed for late B/A-type stars.We take Ṁ = 10 −10 M ⊙ /yr, v ∞ = 1000 km/s, which gives R A ≈ 14 R * .Note that this value is very likely a lower limit to the true Alfvén radius since the value used for the surface magnetic field is actually a lower limit and that used for v ∞ is likely an upper limit to the respective true values (R A increases with increasing surface magnetic field strength, and decreases with increasing v ∞ , see Eq. 4).Thus R A > R K for this star as well, and, like HD 175362, it possesses a centrifugal magnetosphere.
Finally, the star HD 36644 has the least known stellar parameters.Chojnowski et al. (2019) reported the star to have a very strong magnetic field with a magnetic field modulus of 18 kG, which serves as the lower limit to the actual maximum surface magnetic field strength.Shultz et al. (in prep.)estimated the star's bolometric luminosity and effective temperature to be log(L/L ⊙ ) = Left panels: periodograms for the full dataset (top, purple), after pre-whitening with the rotational frequency and the first harmonic (bottom, magenta), and after pre-whitening with all harmonics (bottom, black).The rotational frequency is indicated with a red dash, harmonics with blue dashes.The green dot-dashed curve is a polynomial fit to the fully pre-whitened frequency spectrum, used to evaluate the noise floor; the green solid curve shows 4× the noise floor, the threshold for frequency significance.Right panels: the light curve folded with the rotational frequency (top), and residuals (bottom) after subtraction of the harmonic model determined from frequency analysis (red curve).Different sectors are indicated by colour.
1.3 and log(T eff /K) = 3.95 respectively, which gives a stellar radius of 1.9 R ⊙ .
To estimate the rotational period, we acquired the High-Level Science Product light curves from the Transiting Exoplanet Survey Satellite (TESS) via the Mikulski Archive for Space Telescopes (MAST).HD 36644 was observed during sectors 6 (with 30-minute cadence) and 43, 44, and 45 (with 10-minute cadence).The TESS light curve is shown in Fig. 4.
Inspection of the individual sectors revealed long-term variability, suggesting a rotation period longer than the ∼ 1 month duration of each sector.Determining the period therefore requires combination of the data from multiple sectors.Since the flux in any given sector is normalized by the mean flux in that sector, this leads to systematic offsets in the flux levels.
An iterative approach was taken to determine the period.First, the approximate period was determined using the frequency analysis package period04 (Lenz & Breger 2005).The light curves were then folded with this period, and individually adjusted by adding an ar-bitrary flux offset in order to bring them into alignment.The frequency analysis was then repeated in order to obtain a refined period and a harmonic model for the variation by inclusion of all significant frequencies in the periodogram (where 'significant' is defined as 4× the noise floor).The harmonic model was then used to perform an empirical detrending, by subtracting the harmonic model from each sector, fitting a low-order polynomial to the residuals, and then subtracting the residuals from the original light curve.These de-trended light curves were then used to obtain the final period.
The results are shown in Fig. 4. The light curve displays a double-wave variation, with the strongest peak in the frequency spectrum at 2f rot , together with two additional significant harmonics of f rot .The rotational frequency is f rot = 0.0210638(3) d −1 , corresponding to a rotational period of P rot = 47.4746(6)d.
If this rotational period is correct, the corresponding Kepler radius is 36 R * assuming a mass of 2 M ⊙ .If we assume the other parameters to have the same value as those assumed for HD 40759, the extremely strong sur-face field implies that the Alfvén radius is at least 39 R * .Thus, despite the very long rotation period, it is likely that the star has a thin centrifugal magnetosphere.We emphasize that the lower limit for the Alfvén radius is likely very conservative, first because 18 kG is the lower limit for the surface magnetic field strength, and second because the mass-loss rate of 10 −10 Ṁ yr −1 is probably considerably higher than the actual mass-loss rate of this Ap star, which may easily be orders of magnitude lower.In the future, it will be important to obtain more accurate parameters for the star so as to confirm the presence (or absence) of a centrifugal magnetosphere surrounding the star.

Flare energies
The flux density of HD 175362 at 50 cm is ≈ 0.54 mJy (Shultz et al. 2022), corresponding to an isotropic spectral luminosity (defined as 4πd 2 S, where d is the distance to the star and S is the observed flux density) of ∼ 10 16 ergs/s/Hz (using the distance to the star as 153 pc, Gaia Collaboration et al. 2018).The isotropic spectral luminosity corresponding to the VCSS flux density is ∼ 2 × 10 18 ergs/s/Hz, thus higher by a factor of 100 than the incoherent spectral radio luminosity.
In the case of HD 40759, the flux density of the flare observed from its direction is around 100 mJy, this corresponds to an isotropic spectral luminosity of ≈ 2×10 19 ergs/s/Hz (assuming stellar distance of 430 pc, Gaia Collaboration et al. 2018).Unfortunately, the star's basal radio luminosity is unknown so far, but the inferred spectral luminosity is much higher than that obtained for incoherent radio emission from any hot magnetic star (e.g.Leto et al. 2021).From the VLASS Quick Look image local noise map, the 3σ upper limit to the star's flux density at 3 GHz is 0.6 mJy.
For HD 36644 also, the basal radio flux density is unknown.We obtained a VLASS 3σ upper limit of 0.4 mJy for its 3 GHz flux density.Using the stellar distance as 384 pc (Gaia Collaboration et al. 2018), the isotropic spectral luminosity corresponding to the flare detection (∼ 100 mJy) is ≈ 2 × 10 19 ergs/s/Hz, which is again too high to be of incoherent origin.

Rotational phases corresponding to the flare detection
The only known mechanism for the production of coherent radio emission by hot magnetic stars is ECME, which gives rise to radio pulses at rotational phases around the magnetic nulls (e.g.Trigilio et al. 2000).In this subsection, we examine the stellar rotational phases corresponding to their radio detections.
The HJD corresponding to the VCSS detection of HD 175362 is 2459619.19268.Using the ephemeris of Shultz et al. (2018), we obtain the corresponding rotational phase as 0.16 ± 0.02.This rotational phase is significantly offset from the magnetic nulls at phases 0.38 ± 0.01 and 0.77 ± 0.01, but close to a ⟨B z ⟩ extremum at phase 0 (Fig. 5, top).In the case of HD 40759, the ⟨B z ⟩ curve never passes through zero though it approaches it very closely (Semenko et al. in prep.).Using the ephemeris: HJD = 2458487.0173 + 3.37484 E (Semenko et al. in prep.), we obtain the rotational phase corresponding to our radio detection as 0.84.Interestingly, phase 0.84 is again close to a magnetic field extremum (Fig. 5, bottom).
The VCSS detection of HD 36644 occurred on HJD = 2458028.8900.Using the ephemeris from our TESS analysis in Sect.5.1: HJD = 2459544.8(1)+ 47.4746(6) E, we obtain the rotational phase as 0.07 ± 0.01.Unfortunately, the correspondence to the longitudinal magnetic field cannot be determined due to the current paucity of magnetic data for this star.5.4.Can the flares originate at the stellar magnetospheres?
Under the scenario that the VCSS detections corresponds to flares from the stars, the underlying radio emission is almost certainly coherent in nature.Note that for CU Vir, in addition to the transient enhance-ments observed close to the extrema of the longitudinal magnetic field ⟨B z ⟩, Das & Chandra (2021) also observed a flux density enhancement of a magnitude (relative to the basal flux density of the star) similar to that observed for the stars reported here.This enhancement was observed at a rotational phase where regular pulses produced by ECME were expected.Hence, this was attributed to ECME, and the unusual strength of the enhancement was speculated to be related to an unusually strong CBO event (Das & Chandra 2021).
For the case of a star with an axi-symmetric dipolar magnetic field, ECME pulses ideally appear around the rotational phase corresponding to magnetic nulls.For the two stars for which we could correlate the rotational phases corresponding to the VCSS detections, the rotational phases are greatly offset from the phases of minimum |⟨B z ⟩|, and they rather lie close to the maxima of |⟨B z ⟩|.This situation is similar to the case of the radio flares observed from CU Vir (Das & Chandra 2021), and supports the idea that they all represent a class of coherent transient events different from the periodic ECME phenomenon.
For the adopted distances, the estimated radio luminosities (∼ 10 18 − 10 19 ergs/s/Hz at 340 MHz) are some of the highest ever observed for stellar flares, and rules out certain alternative sources.M-dwarfs are one of the most common sources of stellar flares.From Figure 1  The only sources that produce flares above a luminosity of 10 17 ergs/s/Hz (at 144 MHz) are the RS CVn systems and FK Com stars.So far, none of the stars are known to have any such companions.However, Bodensteiner et al. (2018) reported the detection of an infrared nebula around HD 175362, which they attributed to binary interactions.Thus, it is important to perform detailed follow-up observations of all three stars in the future in order to characterize their radio properties, and also at other wavebands to search for potential companions.
In addition to being possible sources of the radio flares themselves, companions can also play an indirect role in the generation of flares from hot star magnetospheres in a manner similar to how Io triggers auroral radio emission from Jupiter (e.g.Belcher 1987).In this case, when the companion moves through the magnetosphere of the hot star, it may lead to enhancement in the non-thermal electron density and manifest as a radio flare.
While flares of CBO origin will be random in terms of their occurrence times, those triggered by a compan-ion are expected to exhibit a periodicity related to the orbital motion.With a single detection, however, it is not possible to distinguish between different scenarios.Follow-up observations must be performed, both for understanding the temporal properties as well as to detect companion(s).One of the key challenges for the latter will be to resolve objects that are only a few milliarcseconds apart.In that case, one can perform Very Long Baseline Interferometry (VLBI) with telescopes such as the High Sensitivity Array (HSA) that offer angular resolution of ≲ 1 mas.Even if the companion is not sufficiently radio bright during non-flaring phases, its presence might be detectable through positional shift of the hot magnetic star due to its orbital motion.

SUMMARY
In this work, we report radio flares at 340 MHz from the directions of three hot magnetic stars using VCSS data.By characterizing the imaging artifacts observed in VCSS snapshots, we find the probability that these three detections are false associations with artifacts is less than 1%.This low probability motivates us to examine whether the stars can produce radio flares in their magnetospheres.Using the available stellar parameters, we conclude that all three stars are suitable for harboring centrifugal magnetospheres where CBOs can occur.Note that HD 175362 is confirmed to have a centrifugal magnetosphere and HD 40759 is extremely likely to have one.HD 36644 needs to be investigated further to infer if it indeed possesses a centrifugal magnetosphere.We also find that similar to CU Vir (the only other hot magnetic star from which radio flares have been reported) the flares were observed close to a magnetic extremum, for the two stars with sufficient magnetic data.Whether this has implications for the underlying flare production mechanism can only be confirmed by more such detections.
To summarize, although our current data do not allow us to confirm the source of the observed enhancements in the radio light curves, neither does the data suggest the stars cannot be the origin of the radio flares.It will be interesting to see if more flares are detected from hot magnetic stars in the ongoing VCSS epoch 3, or by other radio surveys, since statistical analysis is probably the only way to confirm whether hot magnetic stars flare or not.From these results we conclude VCSS artifacts can be identified as sources without TGSS or RACS catalog matches and that the nearest neighbor distance can be used to determine which population artifacts belong to.We use this method to count the total number of isotropic artifacts with SN R > 5 in each snapshot.
The right panel of Figure 6 shows the isotropic artifact count distribution in snapshots within either the TGSS or RACS footprint.Over 99.7% of snapshots in each epoch are within this footprint.There is a large variation in the number of artifacts per snapshot with distribution peaks at 5 − 10 artifacts.

B. FALSE ASSOCIATION SIMULATIONS
To confirm our statistical calculations we ran simulations generating 50, 000 random catalogs using a model for the observed sky distribution.We assume the distribution is symmetric about the Galactic Plane and plot the stellar density binned in galactic latitude in the left panel of Figure 7.We fit the density with an exponential function:  The middle panel of Figure 7 compares the catalog normalized galactic latitude cumulative density function (CDF) to a numerical integration of the fitted stellar density function.The fit is a good approximation to the data.The right panel of Figure 7 shows the CDF of a uniform distribution is a fair approximation to the catalog distribution in galactic longitude.Mock catalogs are constructed using an inverse transform technique to generate random coordinates for 761 stars following the model distribution.A sample from a random variable uniformly distributed between 0 and 1 is generated.The star's galactic latitude is set to where the CDF equals this number.Another uniform random variate is used to set the galactic longitude assuming stars are uniformly distributed along this axis.Figure 8 compares the distributions of the stellar catalog and a mock catalog.
Figure 9. Fraction of 50,000 simulated catalogs having N match matches to the VCSS catalog.A Poisson distribution with expected value equal to the average number of matches (0.34) is a good fit to the simulations (dashed line).The simulations agree with the calculations for the stellar catalog (solid line) that the probability of three VCSS matches is < 1%.
The number of VCSS matches is calculated for each mock catalog.The fraction of catalogs having a given number of matches is shown by the blue points in Figure 9.The distribution is well described by Poisson statistics with an expected value equal to the average number of matches, λ = 0.34 (dashed line).Only 238 mock catalogs have three VCSS matches, giving a false association probability of 0.5%.This is lower than our calculation for the stellar catalog (solid line) and likely due to assumptions made on the stellar distribution: symmetry about the Galactic Plane and smoothly distributed in longitude (no clumps); as well as uncertainty in our fit to the density function.The simulations agree, however, that the probability of all three VCSS hot magnetic star detections are due to false associations with artifacts is < 1%.

Figure 1 .
Figure 1.Light curves of hot magnetic stars detected in VCSS snapshots.The time axis is aligned with the UT of each detection, indicated and labeled in each plot.Upper limits for non-detections (triangles) correspond to three times the PyBDSF-calculated local noise at the stellar positions corrected for primary beam attenuation .

Figure 2 .
Figure 2. VCSS snapshot cutouts centered on HD 36644, organized chronologically.Pixel intensities have been corrected for primary beam attenuation.Labels indicate the time relative to the detection, the distance of HD 36644 from the snapshot center and the noise (σ) at the stellar position (circled).The color scale ranges 0-6 times the local noise for each image.

Figure 3 .
Figure 3. VCSS snapshot cutouts centered on HD 40759 (top) and HD 175362 (bottom) at the time of detection.Color scale is as in Figure 2.

Figure 4 .
Figure4.TESS photometry of HD 36644.Left panels: periodograms for the full dataset (top, purple), after pre-whitening with the rotational frequency and the first harmonic (bottom, magenta), and after pre-whitening with all harmonics (bottom, black).The rotational frequency is indicated with a red dash, harmonics with blue dashes.The green dot-dashed curve is a polynomial fit to the fully pre-whitened frequency spectrum, used to evaluate the noise floor; the green solid curve shows 4× the noise floor, the threshold for frequency significance.Right panels: the light curve folded with the rotational frequency (top), and residuals (bottom) after subtraction of the harmonic model determined from frequency analysis (red curve).Different sectors are indicated by colour.
of  Vedantham et al. (2022), M-dwarf flares have luminosities below 10 15 ergs/s/Hz at 144 MHz.Thus M-dwarf companions to the stars are unlikely to produce flares such as the ones detected by VCSS.

Figure 6 .
Figure 6.Left: Average source density around isolated point sources in VCSS epoch 1 (dashed) and epoch 2 (solid).Densities above several SNR cuts are shown.Right: Histograms of the number of isotropic artifacts in VCSS snapshots.

n
= n 0 e −k|b| + c (B3)where n is the star density per steradian, b is the galactic latitude in radians, and n 0 , k and c are constants.The fitted values and 1σ uncertainties are (n 0 , k, c) = (191 ± 12, 3.8 ± 0.5, 13 ± 4).The two bins at highest declination were not used in fitting due to large uncertainties from the paucity of sources.

Figure 7 .
Figure 7. Left: Surface density with galactic latitude of the hot magnetic star catalog and the fitted equation (B3).Middle: Normalized galactic latitude CDF for the catalog and fit.Right: Normalized galactic longitude CDF for the catalog and a continuous uniform distribution.

Figure 8 .
Figure 8. Left: Catalog of 761 hot magnetic stars plotted in Galactic coordinates.Right: A randomly generated catalog with the same number of stars and similar distribution.