From downtown to the outskirts: a radio survey of the Orion Nebula Cluster

We present a newly enlarged census of the compact radio population towards the Orion Nebula Cluster (ONC) using high-sensitivity continuum maps (3-10 $\mu$Jy bm$^{-1}$) from a total of $\sim30$ h centimeter-wavelength observations over an area of $\sim$20$'\times20'$ obtained in the C-band (4$-$8 GHz) with the Karl G. Jansky Very Large Array (VLA) in its high-resolution A-configuration. We thus complement our previous deep survey of the innermost areas of the ONC, now covering the field of view of the Chandra Orion Ultra-deep Project (COUP). Our catalog contains 521 compact radio sources of which 198 are new detections. Overall, we find that 17% of the (mostly stellar) COUP sources have radio counterparts, while 53% of the radio sources have COUP counterparts. Most notably, the radio detection fraction of X-ray sources is higher in the inner cluster and almost constant for $r>3'$ (0.36 pc) from $\theta^1$ Ori C suggesting a correlation between the radio emission mechanism of these sources and their distance from the most massive stars at the center of the cluster, for example due to increased photoionisation of circumstellar disks. The combination with our previous observations four years prior lead to the discovery of fast proper motions of up to $\sim$373 km s$^{-1}$ from faint radio sources associated with ejecta of the OMC1 explosion. Finally, we search for strong radio variability. We found changes in flux density by a factor of $\lesssim$5 within our observations and a few sources with changes by a factor $>$10 on long timescales of a few years.


INTRODUCTION
The advent of wideband centimeter-wavelength observing capabilities has enabled a new era of stellar radio astronomy, including observations of radio counterparts of young stellar objects (YSOs). Generally, these show evidence for both thermal free-free emission from ionized material and also nonthermal (gyro-)synchrotron emission from magnetospheric activity (Dulk 1985;Güdel 2002). A few years ago, in Forbrich et al. (2016), we used the National Radio Astronomy Observatory's upgraded Very Large Array (VLA) to study radio counterparts of YSOs in the Orion Nebula Cluster (ONC), the richest nearby young stellar cluster. This study, based on pointing the VLA at the heart of the ONC for about 30 hours, increased the number of known radio sources in the ONC by a factor of ∼7 with the detection of 556 compact sources. Other studies, using the VLA toward the ONC have also focused in the inner cluster covering similar areas (within ∼ 6 × 6 ) usually reaching rms noise levels above 30 Jy bm −1 limiting the number of detections ★ E-mail: j.i.vargas-gonzalez@herts.ac.uk from a few tens (prior to the VLA upgrade Churchwell et al. 1987;Garay et al. 1987;Felli et al. 1993;Zapata et al. 2004) to up to 175 sources (Sheehan et al. 2016). In a larger surveyed area, Kounkel et al. (2014) obtained a shallow map of approximately 1. • 6 × 0. • 4 around the ONC reporting a total of 165 sources and typical rms noise levels of 60 Jy bm −1 .
Additionally, multi-epoch VLA data with its high angular resolution (0. 1) and astrometric capabilites have been used to constrain the kinematics of the ONC (Gómez et al. 2005;Kounkel et al. 2014;Dzib et al. 2017). The main focus has been on the improvement of proper motion (PM) estimates of the main stellar radio sources in the Orion Becklin-Neugebauer/Kleinmann-Low region (BN/KL) in the inner ONC (Becklin & Neugebauer 1967;Kleinmann & Low 1967;Gómez et al. 2008;Zapata et al. 2009;Rodríguez et al. 2017Rodríguez et al. , 2020. The PM of these stars supports the scenario where a multiple stellar system experienced a close dynamical interaction resulting in two massive stars (BN source and source I) and at least one lowmass star being ejected at a few tens of km s −1 triggering a powerful outflow emerging from the OMC1 cloud core prominently seen at near-infrared wavelengths (Bally et al. 2015). The improved sen-sitivity of the VLA now enables the detection also of non-stellar emission like that from jets and outflows (Forbrich et al. 2016;Bally et al. 2020), enabling astrometric studies as presented here.
While the deep radio catalog presented in Forbrich et al. (2016) considerably improved the census of compact radio sources it left open the question of the wider radio population in the ONC, and its interplay with the well-characterized X-ray and infrared populations. The Chandra Orion Ultra-deep Project (COUP; Getman et al. 2005b) covers a larger area (∼ 17 × 17 ) than our single, deep VLA pointing around the same reference center, and we thus conducted a wider survey for radio sources in the ONC which is presented in this work. Six additional pointings surrounding the deep central pointing were obtained at an unprecedented sensitivity in this area and we also repeated the central pointing for comparison (see Figure  1).
The new observations discussed here allow us to obtain the largest census to date of radio counterparts to YSOs anywhere, which we can place into the rich multi-wavelength context of the ONC. We additionally make use of radio astrometry in a comparison of the two central pointings, separated by only ∼4.15 years, to study fast proper motions in the ONC, which is mainly of interest for nonstellar emission, which otherwise are more difficult to measure while providing valuable additional information for source identification. Finally, we use the excellent sensitivity even on short timescales to continue our study of YSO radio variability, motivated by the findings of extreme variability on short timescales (factor >138 in less than two days and a factor of 10 in less than 30 minutes) in our previous deep ONC pointing (Forbrich et al. 2017).

OBSERVATIONS AND DATA REDUCTION
The radio data were obtained between October and November 2016 using the NRAO 1 VLA (project code: 16B-268). Figure 1 shows the observational setup which consists of a central pointing at ( , ) 2000 = (5 h 35 m 14 · 5, −5 • 22 m 30 · 6) and six adjacent pointings observed for about 4 hrs each with the most extended A-configuration array in C-band. The phase centers and dates of the observations are listed in Table 1. The half-power beamwidth (HPBW) for the low-frequency limit (4288 MHz) is ∼ 10. 5 covering a total area of ∼ 20 × 20 arcmin within the collective low-frequency HPBW of all the pointings. The receivers were in full polarization mode with two basebands of 1 GHz each centred at 4.8 and 7.3 GHz with a total of 16 spectral windows (8 per baseband) divided into 64 channels of 2 MHz width each. The primary flux density calibrator for all the pointings was 3C48 and the phase/gain calibrator was J0541-0541 observed every 5−6 min to ensure phase stability as in our earlier observations. The light-blue circles in Figure 1 represent the HPBW at the low (4288 MHz) and high (7847 MHz) frequency ends of the bandwidth in dashed and continuous lines (∼ 10. 5 and ∼ 5. 8, respectively).
The reduction of the data was performed using the VLA Calibration Pipeline using the CASA 2 software (release 5.4.1). All pointings except for pointing 5 were reduced with the automatic processing of the pipeline. Pointing 5 required additional manual selection of faulty data to be excluded in a small portion of the observation which includes 4 science scans and 1 calibrator scan 1 National Radio Astronomy Observatory is a facility of the National Sci- (equivalent to a 5 min interval). No time or spectral averaging was applied to any of the different pointings.
The calibrated data were imaged with the TCLEAN task in CASA. All pointings were imaged to have a size of 8192 × 8192 pixels with a pixel size of 0. 1 as a compromise between the pixel coverage of the synthesized beam and the final size of the image in order to cover the largest half-power beamwidth of ∼10. 5. We used the Stokes plane and spectral definition mode 'mfs' (Multi-Frequency Synthesis) that combines the data from all the selected spectral channels into a single continuum image. The Hogbom deconvolution algorithm and a Briggs weighting method with a robustness parameter of 0.5 were used. For the central pointing, an additional set of images was created with similar setup but using the 'mtmfs' algorithm (Multi-term Multi-Frequency Synthesis; Rau & Cornwell 2011) with 'nterms=2', generating Taylor-coefficient images corresponding to a continuum intensity and spectral index map. This imaging method used for the central pointing enables a direct comparison with the deep radio observations described in Forbrich et al. (2016). Spatial filtering of the visibility data was applied using baselines of the ( , ) range longer than 100 k (∼6 km) to reduce the impact of extended nebular emission on the extraction of point-like sources by filtering out structures greater than ∼ 2 .
The main parameters of the resulting images for each pointing are listed in Table 1, including the synthesized beam sizes and noise levels. The synthesized beam sizes are typically around 0. 3 with only a slightly larger major axis of ∼ 0. 5 for pointing 5. The rms noise level, on the other hand, varies considerably throughout the different pointings reaching the highest value in the central pointing due to the complex structure of the inner part of the cluster. All the adjacent pointings have rms noise levels around 3 − 5 Jy bm −1 reaching the lowest value in pointing 1 in the north-west. The deep observation of the central pointing presented in Forbrich et al. (2016) has a nominal rms noise of 3 Jy bm −1 being the most sensitive observations of the inner ONC to date at these frequencies.
We reached similar rms noise in the outer pointings where most of the crowded area and complex structures lie towards the edges on these images. However, they still represent the most sensitive observations of the ONC to date while similar studies in this region have reported rms noise levels in the range of 25−80 Jy bm −1 (Zapata et al. 2004;Kounkel et al. 2014;Sheehan et al. 2016). To account for the radial decrease in sensitivity caused by the wideband primary beam response we have applied a primary beam correction factor to each source in our catalog after imaging the data following the method described in Forbrich et al. (2016). This primary beam correction factor is a function of distance to the phase center described by a polynomial.
For direct comparison with our earlier results and in order to determine proper motions, while widefield imaging within CASA is preliminary, we have used the standard gridder for our images.
Our imaging experiments show that presently available widefield imaging (w-projection method through gridder='wproject') results in standard-gridder positions in the outer beam (r> 6. 4) that can be slightly off by up to 0. 36 from the corresponding wproject positions, i.e., up to the size of the synthesized beam, however there is only one source at such a large distance in our catalog presented in section §3.1. At distances < 3 from the phase center (∼80% of the sources in our catalog) the offsets found between these two gridders are negligible at 0. 04 ( 10% the size of the synthesized beam). For even greater distances, the offset rises nonlinearly to about half the synthesized beam size (0. 2) at = 5. 2, encompassing 97% of our sources. As discussed below, in our paper this has only a minor  Forbrich et al. (2016), is surrounded by six additional pointings (listed in Table 1) with the same spectral setup in C-band. The light-blue circles indicate the HPBW of each pointing at the low (4288 MHz) and high (7847 MHz) frequency ends of the bandwidth in dashed and continuous lines (∼ 10. 5 and ∼ 5. 8, respectively). Red symbols show the radio sources detected in this work, green symbols indicate radio sources detected in the deep survey (Forbrich et al. 2016), with yellow symbols additionally marking the positions of 1 Ori C and the BN object for reference (lower left and upper right, respectively). Blue symbols indicate the positions of X-ray sources from the COUP survey (Getman et al. 2005a). The background image is a HST -band image (ACS/WFC) of the Orion Nebula (Credit: NASA, ESA, M. Robberto, and the Hubble Space Telescope Orion Treasury Project Team). Note: The light-blue number labels in Figure 1 correspond to the pointing numbers in this Table. effect in catalog cross-matching for a limited number of sources in our catalog. The imaging applied in both VLA epochs is identical, except that for the new epoch we have not applied any channel averaging, which is different from the approach in Forbrich et al. (2016). In order to quantify any impact from using these two different approaches we re-imaged the central pointing of the new observations with the same channel averaging used before. The positions in both images are compatible within the uncertainties where even at a large distance from the phase center (r=6 ) the effect corresponds to a shift of less than 0.5 (∼22 mas). The impact on flux densities is equally negligible.

Source Detection and Distribution
The high resolution provided by the A-configuration array together with the additional spatial filtering of the visibility data applied in the imaging process largely mitigated the difficulties of disentangling the compact radio emission from the range of emission size scales over the area due to the complex structure of the Orion Nebula. However, as reported by Forbrich et al. (2016), the use of automated methods for point source extraction leads to a high fraction of spurious detection (up to 50%) on VLA images with the same observational setup towards the ONC. The source detection was thus performed by visual inspection of each individual pointing using the images generated in the multi-frequency synthesis (mfs) spectral mode. The initial source positions in this process were estimated with DAOFIND task in IRAF. This task computes the positions by estimating the point spread function of the source using an elliptical Gaussian approximation within a given area defined by a box around the source. This first list of detections per pointing corresponds to the input for a source extraction script which uses the IMFIT task in CASA to obtain the final positions, flux densities and additional statistical parameters per source. For each input source we used different box sizes for the fitting ranging from 8 to 30 pixels, applying different offsets from the center of the source to avoid contaminant emission coming from nearby sources. The different outputs for a given source from different box configurations were compared to finally select the measurement with lower uncertainties and measured positions closer to the actual peak pixel, thereby minimising the impact of nearby sources and of complex nebular emission. We also enforced a minimum box size of 10×10 pixels for it to contain at least about ten synthesized beam areas for statistics. To select reliable detections we used a S/N>5 as our detection limit. Compared with a criterion of S/N>3, our adopted limit conservatively accounts for non-Gaussian noise in the inner cluster, which leads to many more sources that would need to be rejected (e.g., compact nebular emission). A 5 cutoff is also consistent with the most complete radio surveys to date in the ONC field with detection limits between 4.5 − 6 (Kounkel et al. 2014;Sheehan et al. 2016;Forbrich et al. 2016).
The final number of sources detected per pointing is listed in column (2) in Table 2 and marked in red in Figure 1. Column (3) in the table indicates the number of sources within the HPBW (at the high frequency end ∼ 5. 8) for each pointing. The maximum number of detections is found in the central pointing with 272 detections (which corresponds to the densest region of the ONC) and the lowest detection count is found in pointing 5.
The total number of detections over the whole area covered by the 7 pointings is 521, without duplicates that are detected in more than one pointing. Since most of the duplicates are isolated sources, those were easily found with a search radius of 1", followed by visual inspection.
Sources detected in several pointings are usually found at different distances from the phase centers. As noted above, these different detections show minor position shifts. In our catalog we thus report positions from the closest detection to the phase center, where this corresponding distance is also reported. Partly due to the spatial filtering inherent in our experiment design, focusing on the detection of nonthermal emission with aggressive spatial filtering, most of the sources are unresolved, with only 2% showing a size 3 times the area of the synthesized beam, and only 8% of the sources larger than twice the synthesized beam. The nominal mean ratio between the integrated to peak flux density in our catalog is 1.2 ± 0.5. Our catalog thus only lists peak flux densities.
Our catalog of 521 sources is listed in Table 3 and indicated by red symbols in Figure 1. Columns (1) and (2) show the positions in and with their corresponding uncertainties obtained from the fit (IMFIT). Column (3) indicates the source identification number in this catalog. Columns (4) and (5) indicate the peak flux density (corrected by the primary beam response) and source fitting parameters (major and minor axes, and position angle). Columns (6) and (7) indicate previous designation in the COUP and/or VISION surveys. Column (8) indicates the distance to the stellar system 1 Ori C and column (9) indicates the distance to the closest phase center where the reported positions come from. The position uncertainties reported in this catalog are those given by IMFIT without the addition of minor systematic errors (see above) 3 . These positional uncertainties have median values of 14 and 15 mas in R.A. and decl., respectively. An absolute astrometric accuracy for similar VLA observations (identical to our central pointing) was reported in Forbrich et al. (2016) using five individual epochs resulting in an overall absolute astrometric accuracy of 20-30 mas. An absolute uncertainty in peak flux densities of 5% has been estimated based on systematic variability using a non-variable test-case (source BN; see section §3.4). This uncertainty has been added in quadrature with the uncertainties from the 2D-Gaussian fit (IMFIT) already corrected by the primary beam response.
In this new census of compact radio sources towards the ONC we report 198 new sources not previously reported at these frequencies. This is based on a comparison against the most complete catalogs of compact radio sources made of the ONC and presented in  As defined by imfit in CASA. Distance to closest phase center. Note that the reported positions are from the standard gridder and may be slightly off in the outermost beam areas (see text). The full catalog is available as supplementary material.

Multi-wavelength Populations
Our primary goal after revealing the compact radio population in the ONC is the search for radio counterparts to YSOs. They can be detected through both thermal and nonthermal radio emission. While nonthermal emission originates at the smallest scales in the stellar coronae of low-mass young stars, thermal radio emission, in contrast, occurs in a range of larger scales in stellar disks or outflows when it is associated to a stellar source. However, compact sources at these frequencies can also be detected as thermal emission from ionized material not strictly related to a stellar source (Garay et al. 1987;Churchwell et al. 1987;Sheehan et al. 2016;Forbrich et al. 2016). An additional source of contamination comes from extragalactic background sources, even though these have been previously found to be minimal. In the so far deepest catalog for the inner ONC (Forbrich et al. 2016) it was estimated that ∼97% of the compact radio sources within < 1. 6 are related to the cluster. Here we cover a considerably wider area and while the central region is less prone to present background detections due to the elevated rms noise caused by the nebula itself, we expect to find faint background radio sources from bright background galaxies in the outer areas, given the high sensitivity. A discussion on the expected number of extragalactic radio sources in the outer areas is presented later on in this section, when discussing Figure 5.
In order to approach the nature of the compact radio population we have compared our catalog against the most complete X-ray and near-infrared (NIR) datasets available for the ONC. The COUP Xray catalog reported in Getman et al. (2005b) represents the best reference for YSOs in the ONC. Here, the X-ray emission of young stars is associated to thermal emission from hot plasma in coronaltype activity (e.g., Feigelson & Montmerle 1999). The observations presented in this work were designed to cover the COUP survey area of ∼20 ×20 (blue markers in Figure 1 represent the COUP catalog). The complementary NIR dataset used here is the VISION catalog reported in Meingast et al. (2016), a survey of the entire Orion A molecular cloud in bands. The NIR band is an additional tracer of the young stellar population towards the ONC, since, in general, such X-ray and NIR data are similarly effective at detecting embedded sources (e.g., Ryter 1996). Any extragalactic background contamination is unlikely where sources are superimposed onto the high extinctions levels of the ONC but at the same time the NIR also picks up foreground sources. Contrary to the X-ray survey, the NIR catalog is severely affected by bright extended emission of the Orion Nebula in the innermost areas, where it is thus less complete to YSOs.
The angular resolution in COUP and VISION is < 0. 5, which is comparable to the typical beam size in our observations. Moreover, the reported positional errors in the COUP survey have a mean value of 0. 1, with maximum errors of 0. 59, although, 99% of the sources have positional errors 0. 5. On the other hand, the reported astrometric accuracy in the VISION catalogue is ∼70 mas. Based on these considerations and to additionally account for the smallest nearest-neighbour distances, we have therefore considered a conservative radius of 0. 5 for the search of counterparts in these surveys. The robustness of this approach is demonstrated when considering that if we increase this search radius to 1 the number of correlations with the COUP and VISION surveys would increase by only 5% and 4%, respectively, including unrelated nearest neighbours. Within the area covered in our observations we report 521 radio sources, while the COUP catalog reports 1616 sources and the VISION catalog 3558 sources. A total of 275 sources in our catalog (∼53%) have X-ray counterparts which corresponds to only 17% of the COUP catalog. On the other hand, 290 NIR sources (or ∼56%) have counterparts in our radio catalog, which corresponds to only ∼8% of the NIR catalog within the same area. Figure 2 shows the surface number density of the three wavelength populations as a function of projected distance from 1 Ori C, a young and massive stellar system of ∼50 in the Trapezium cluster at the center of the ONC (Kraus et al. 2007). Our obser- vational setup is actually centered ∼1 northwest from 1 Ori C towards the BN/KL region to optimize the sensitivity in this complex area while still matching the COUP field of view. Centering our reference on 1 Ori C allows us to address any impact that the most massive stars in the center of the cluster have on the radio emission in a multiwavelength context. The overall distributions for the three bands show the densest region in the central cluster ( < 1 ) with a remarkable correspondence in their total number of sources equivalent to an average surface density of (3.7 ± 0.2) × 10 3 sources pc −2 . Despite of this similarity, they do not represent the same population and, indeed, at this inner bin their actual correlation, indicated by the orange distribution (population with detections in the three bands), is only one third of their total number of sources. The X-ray and NIR distributions exhibit almost the same profile until = 2 . The X-ray distribution continuously decreases revealing, in part, the structure of the cluster but also showing the sensitivity limitations in the outer bins. Similarly, the radio distribution continuously decreases, although at a faster rate. Contrary to the deep radio catalog distribution (dashed black line), our radio catalog (continuous red line) does not decrease solely due to a sensitivity effect, even though the sensitivity is not constant. The angular distance between the phase centers of adjacent pointings is ∼ 5 where the lowest sensitivity occurs at 2. 5 between each pointing. At this distance, the wideband primary beam correction indicates a flux density correction by 40%, which therefore technically constitutes the maximum variation in sensitivity due to the spacing of our individual pointings. However, we have shown in Forbrich et al. (2016) that at this angular distance from the phase center of the central pointing the image noise is still dominated by the Orion Nebula itself, and the impact on our analysis is thus limited. Under this considerations, the radio distribution seems to reveal the intrinsic structure of the cluster. At larger radii ( > 4 ) the NIR distribution remains almost constant largely due to the presence of foreground infrared sources in addition to the young stellar population further away from the center of the cluster. Although the COUP distribution may also include contamination, in this case from extragalactic candidates, this is just a small fraction (∼160) and most of the Xray sources are likely members of the ONC (Getman et al. 2005b), besides, extragalactic X-ray sources will be affected by foreground extinction in the cloud and thus are unevenly distributed.
In order to assess the correlation between the three main distributions shown in Figure 2 (radio/X-ray/NIR) we quantified the different population fractions as a function of projected distance from 1 Ori C as shown in Figure 3. Since we have shifted our reference center ∼ 1 southeast from the actual center of the observations, we have therefore moved the reference frame for this analysis towards the south-east boundary of our radio survey (and COUP given the almost identical coverage). Since the maximum radial coverage in our radio catalogue is ∼ 10. 4, we are now reaching this boundary at ∼ 9. 4 south-east from 1 Ori C, setting a first limit for our analysis at this maximum radius. However, an additional constraint is set to avoid any sensitivity bias for the detection of sources in the boundaries of our radio catalogue (which becomes sensitive to the brightest sources due to the primary beam correction). We have set a radial limit for our analysis to not exceed the midpoint between the HPBW at the low-and high-frequency ends of the bandwidth (dashed and continuous line circles, respectively, in Figure 1). This midpoint is reached at 8 in the south-east direction, thus we are going to restrict our following analysis to this radius. The left panel in Figure 3 indicates the fraction between the total radio distribution shown in Figure 2 over the total X-ray and NIR distributions (indicated by blue and green histograms, respectively), without matching individual counterparts. In the central cluster these fractions indicate almost equal numbers of X-ray and NIR sources with radio source numbers lower by ∼ 10−30%. These fractions then decrease as expected from their individual distributions in Figure 2. Here the radio to NIR fraction (green) decreases at a slightly faster rate due to the more uniform spatial distribution of NIR sources in the field (including foreground sources) compared to the continuously decreasing surface density distribution of radio sources as we go further away from 1 Ori C. Interestingly, for projected distances of > 3 , unlike the continuous decrease in the radio to NIR fractions, the radio to X-ray fractions (blue) settles at around 19 ± 1% which might indicate how these two populations are tracing a similar structure of the cluster within this radial range.
The middle panel in Figure 3 shows the X-ray and NIR detection fraction of radio sources (i.e., X-ray or NIR counterparts to our radio catalog) indicated by blue and green histograms, respectively. The overall distributions do not show a significant trend and they are, in effect, compatible with a constant distribution within the errors. Evidently, at larger radii the lower number of radio sources introduces larger uncertainties. A similar comparison is shown in the right panel, but indicating the radio detection fractions of X-ray and NIR sources (i.e., number of X-ray or NIR sources with radio counterparts) in blue and green, respectively. Here, the blue distribution represents our best estimate for the radio detection fraction of YSOs in the ONC and any radial trend with respect to the brightest Trapezium star at the center of the cluster may well provide important information on the underlying emission mechanism dominating at different distances. If YSOs have similar geometries, we would not expect to find any trend here, and the radio to X-ray ratio would instead simply reflect basic YSO properties. However, there is a very clear radial trend suggesting that X-ray sources towards the central part of the cluster are more likely to have a radio counterpart than X-ray sources in the outer areas. This points to a potential impact . Detection fractions between the three different populations radio/X-ray/NIR as a function of distance to 1 Ori C. The left panel shows the fraction of radio sources over the X-ray (blue) and NIR (green) populations. Central panel shows the X-ray (blue) and NIR (green) detection fraction of radio sources. The right panel shows the radio detection fraction of X-ray sources (blue) and NIR sources (green). The 1 error bars were derived from counting statistics (Poisson errors).
of the Trapezium on YSO properties in the ONC. One possibility is that circumstellar disks of YSOs (with a central X-ray source) are externally photoionized by the influence of the Trapezium stars (O'dell et al. 1993;Henney & Arthur 1998;Concha-Ramírez et al. 2020 and references therein), leading to the detection of ionized material as thermal free-free radio emission. In this case, as we go further away from the Trapezium, the detection fraction (blue distribution) may become dominated by nonthermal radio emission intrinsic to YSOs. The detection fraction is as high as 40 ± 6% in the inner bin ( < 0.12 pc) and then decreases down to 17 ± 3% in the third bin ( < 0.36 pc). For distances > 3 (0.36 pc) the radio detection fraction of X-ray sources fluctuates around 11 ± 2% which might represent a baseline for the detection of nonthermal radio emission from the X-ray emitting YSO population without the influence of nearby young massive stars that would otherwise lead to a constant distribution. In contrast, the radio detection fraction of NIR sources (green) decreases continuously with radial distance, but this is expected since it reflects the fact that the NIR catalog includes not only cluster members but also sources that are not related to the ONC.
Following the interpretation of the high radio detection fraction of X-ray sources closer to the Trapezium, more likely due to the detection of externally photoionized circumstellar disks, it is expected that the flux distribution of these sources (thermal freefree emission component) decreases with distance from the ionizing source. Figure 4 shows the peak flux density distribution as a function of projected distance to 1 Ori C clearly showing higher radio fluxes closer to the Trapezium with a decreasing trend more evident for 2. 5. Possible deviations from this trend are, first, the discrepancy between the projected and actual physical separation to the ionizing source where the actual physical separations could be considerably larger. Second, a fraction of the sources could be intrinsically bright nonthermal radio emitters.
An interesting question prompted by these results concerns the nature of the radio population without X-ray counterparts. On the one hand, towards the center of the cluster most of these sources likely are thermal radio sources including externally photoionized disks, stellar jets and outflows, and also compact emission from the . Peak flux density distribution as a function of projected distance to 1 Ori C. Blue symbols indicate sources with X-ray counterparts and grey symbols indicate the remaining sources. Red symbols show the median of the blue distribution per 0. 25 bins (for < 2 ) and 0. 5 bins (for > 2 ). The median error bars represent the 25th and 75th percentiles of the data at each bin.
nebula, but towards the outer areas where we do not find a significant influence of the Trapezium stars we would expect a much lower number of these sources. As noted above, while the central area is less prone to background contamination already due to presence of the nebula, this contamination increases in the outer areas. Additionally, there may be a low number of sources not identified as YSOs if these were not X-ray active at the time of the observations. Indeed, it has been found that X-ray sources with radio counterparts are predominantly sources with high X-ray luminosities and there are even a few cases where nonthermal radio YSOs have no X-ray counterpart at all (Forbrich & Wolk 2013;Rivilla et al. 2015;Forbrich et al. 2016). The top panel in Figure 5 shows the radial distribution of the total radio population, as well as the radio population with and without X-ray counterparts. The radio population without X-ray counterparts peaks in the central area, implying that there is a correlation with the inner cluster and these sources are not only background galaxies that would lead to a constant distribution as we see in the outer area (> 3 ). In this outer area the distribution of radio sources without X-ray counterparts is evenly spread over the field as indicated by black symbols in the bottom panel of Figure  5. In the same bottom panel we included in red symbols the radio sources with X-ray counterparts but flagged as extragalactic candidates in the COUP catalog. We found only 11 of these sources in our catalog (∼2%) out of 159 extragalactic source candidates in the COUP catalog, which are evenly spread throughout the whole field (Getman et al. 2005a). In order to assess the possibility to detect extragalactic sources we made use of the work by Fomalont et al. (1991) to estimate the expected number of background sources in a given area above a given flux density threshold. We considered an area between 3 < < 9 excluding the central cluster and where the distribution of radio sources without X-ray counterparts remains nearly constant around 11±3 sources per bin, totaling 65±8 sources. Following the 5 detection threshold used in our catalog and an average rms level of 4 Jy bm −1 (see Table 1) we find an expected number of background sources of ∼ 155 ± 22 and therefore extragalactic sources are a reasonable explanation for at least some of the radio sources without X-ray counterparts in the outer areas. In order to further explore the sample of 65 radio sources between 3 < < 9 we have looked into their correlations with NIR and optical wavelengths. 11 sources have NIR counterparts in VISION, 9 sources have optical counterparts in Gaia EDR3 (Gaia Collaboration et al. 2016;Gaia Collaboration et al. 2021), and just 4 sources have counterparts in both bands, however the nature of those remains unknown. On the other hand, 8 sources are associated with the OMC1 outflow discussed later in section §3.3, and 4 sources are associated with proplyds found in optical HST data (Ricci et al. 2008), leaving us with 55 unidentified sources in this sample and thus potentially extragalactic candidates.
We have focused our analysis on azimuthally averaged radial trends rather than individual pointings. This is justified by the fact that we find only limited evidence for significant differences between the individual pointings. In order to quantify any substantial differences between them we compared their relative multiwavelength distributions using KS tests. We found that the outer pointings have similar distributions and that the only large difference comes from comparing the central pointing with any of the outer ones (see section §B for details on the KS statistics).
Another interesting aspect is the relative fraction of radio sources that have neither X-rays nor NIR counterparts which is significantly lower in the outer fields. This population most likely represents compact radio emission of non-stellar origin. In section §3.3, we identify one category of such non-stellar radio emission by employing proper motion measurements, identifying sources associated with the OMC1 outflow.

Proper Motions
In our comparison between the observations presented here and the deep catalog reported in Forbrich et al. (2016) we were able to implicitly measure the proper motion of 253 sources. This comparison only involves the central pointing in epoch 2016, these are identical observations towards exactly the same phase center and phase-referenced against the same distant quasar such that this comparison basically represents absolute astrometry. However, we are only covering a time baseline of 4.19 years. The mean position errors in our catalog are ∼20 mas and the estimated astrometric accuracy from the five epochs in Forbrich et al. (2016) is 20-30 mas. In order to detect proper motions above these limits a single source would require a minimum velocity of ∼22 km s −1 , which is unusual for a stellar source and therefore we did not expect to find significant motions in our sample. Stars with velocities >10 km s −1 are considered peculiar sources and above 30 km s −1 are already in the runaway regime (Farias et al. 2020;Schoettler et al. 2020). Stellar proper motions in the ONC are typically 2 mas yr −1 (4 km s −1 ) Kim et al. 2019) and only a handful number of sources have fast proper motions with 30 mas yr −1 (60 km s −1 ) (Gómez et al. 2005(Gómez et al. , 2008Rodríguez et al. 2020) most of them associated to the ejected stars in the BN/KL region (which will be discussed later in this section).
In spite of the above considerations, we surprisingly discover high PMs in our sample with projected velocities of up to ∼373 km s −1 for a distance of ∼ 400 pc to the ONC (Großschedl et al. 2018;Kuhn et al. 2019). This finding provides an additional and important tool for source identification. Figure 6 shows the absolute PM diagram in the cos( ) × plane for the full sample. Sources with motions above 5 in at least one direction or are colour coded by their transverse velocity ( ). It is important to note here that these uncertainties were derived by error propagation from the fitting errors (positional uncertainties) thus these may be underestimated. 48 out of 253 sources fulfil this criterion. The remaining sources are marked in grey. With this first cutoff sources with small uncertainties but still with insignificant PM remain (even with motions above 5 ). An additional cutoff is then applied taking into account the PM dispersion to estimate an intrinsic motion noise in the PM diagram in Figure 6. This additional cutoff is described below.
An important caveat to consider in this sample is the fact that not all of the sources represent a stellar counterpart (see §3.2) and therefore an unrealistic dispersion might be injected into the central part of the diagram caused by apparent motions of extended sources that can be seen as compact radio emission in form of thermal bremsstrahlung coming from ionized material around young stars, a peak emission in top of a bow shock, for instance. This apparent motion can be caused by fading and brightening events of the same feature in between the two epochs. On the other hand, given the short time separations between the two epochs we are only able to reveal significant motions above a few tens of mas yr −1 . However, there is a continuum range of motions spanning from −70 to 60 mas yr −1 in and from −65 to 125 mas yr −1 in .
In order to determine a realistic PM dispersion we constrained a Gaussian fit for each axis in Figure 6 to a sample of sources that (1) are compact, (2) do not lie on a complex or extended emission and (3) do not present a clear motion which is considered as sources with position measurements compatible within 3 between the two epochs. This selection gives us a dispersion that represents the intrinsic motion noise in the PM diagram due to position uncertainties. This distribution is shown in dark-blue histograms in Figure 6 for each axis, together with their Gaussian fit also in dark-blue colour. The full sample distribution is showed in orange histograms. The shaded ellipse represent 5 times the dispersion of the Gaussian fit for each axis in order to exclude ambiguous motions within this area and focus the discussion in the most significant motions after the previous cutoff criteria.
The final sample of significant motions is listed in Table A1. It consists of 22 sources with transverse velocities in the range of ∼95−373 km s −1 with most of them lying along the "finger" shaped features towards the north and northwest part of the BN/KL outflow from the Orion OMC1 cloud core (Zapata et al. 2009;Bally et al. 2015) and also known as "Orion Fingers" (Taylor et al. 1984;Allen & Burton 1993). Figure 7 shows a high-resolution NIR image towards the OMC1 outflow adapted from Bally et al. (2015) with the broad-band filter in red, the 2.12 m H 2 and 1.64 m [Fe ] narrow-bands in green and blue, respectively. The OMC1 outflow has an explosive morphology of molecular material consequence of the dynamical encounter of stars which were ejected with speeds of a few tens of km/s more than 500 years ago (Bally & Zinnecker 2005;Rodríguez et al. 2005). Most of the material seen in near-IR image in Figure 7 seems to have been triggered by this event. This finding proves the sensitivity of these VLA observations to non-stellar radio emission and the utility of proper motions as an additional piece of information in counterpart identification.
The sample of fastest proper motions are indicated in yellow symbols in Figure 7. An example of two of the fastest features detected are shown in Figure 8. These are the sources 153 and 155 and are the northernmost sources seen in Figure 7, with transverse velocities of = 350 − 373 km s −1 and uncertainties on the order of 13 km s −1 . The left panel of Figure 8 shows a close-up of Figure  Additional examples of non-stellar compact radio emission are those associated with the jet-like features around the BN object and shown in Figure 9. The radio continuum map corresponds to the 2012 VLA observations from Forbrich et al. (2016) while the white contours in the close-up around BN and the southwest features correspond to the 2016 VLA observations. These jet-like features are labeled following the discussion in Bally et al. (2020) as E2 and E1 east from BN, and SW1, SW2 and SW3 southwest from BN. The spatial distribution and proper motions of these sources can be interpreted, as further discussed in Bally et al. (2020), to be tracing a low-velocity outflow from BN, although it has been suggested that it may contain a low-luminosity protostar or to be tracing ejected clumps produced by the OMC1 explosion Bally et al. 2020). Both sources SW1 and SW3 are moving away from BN towards the west while E2 is moving away from BN towards the east. The absolute PM of these sources are listed in Table 4 also indicating their equivalent transverse velocities for an adopted distance of 400 pc. The source SW1 (Zapata 11;Zapata et al. 2004) was originally reported by Menten & Reid (1995) and its PM have been well constrained on longer timescales using up to 9 VLA epochs ranging from 1985 to 2018 Rodríguez et al. 2020) reporting PM consistent with our measurements. Additionally, Forbrich et al. (2016) reported spectral index measurement for this source with negative value of −0.3±0.09 suggesting nonthermal emission, which has been actually found for YSO jets (see Section 6 in Anglada et al. 2018 and reference therein). Source E1 is not included neither in Forbrich et al. (2016), since it appears as an extended emission in the epoch 2012, nor in the catalog presented in this work since it is too faint in epoch 2016, however, it is still possible to extract its radio properties in both epochs for the purpose of the discussion in Bally et al. (2020), although with comparatively high positional uncertainties (55-90 mas in and 75-120 mas in ). Interestingly, source E1 and SW2 present motions moving back towards BN, however this may be an apparent motion due to intensity variations of different parts in these features, which also leads to a very high PM for E1. The PM measurements for the sources BN, SW1 (Zapata 11) and E2 (IRc23) are consistent with those reported in Rodríguez et al. (2020).   This work (see Table 3). Source designation in Dzib et al. (2017). Source properties does not meet the criteria for the catalogs defined by neither Forbrich et al. (2016) (it's a extended emission) nor in this work (it has SNR<5).

Peculiar motion of source 117: COUP 510
In the sample of high proper motions there is an intriguing stellar source with apparent fast motion of ( ) = −67.0 ± 1.7 mas yr −1 and = 55.7±2.2 mas yr −1 equivalent to an unusual stellar transverse velocity of = 165 km s −1 . It has an X-ray counterpart in the COUP survey (COUP 510) presenting high extinction with hydrogen column density log( ) = 23.54 ± 0.06 cm −2 (Getman et al. 2005b) and no detection in the VISION Ks-band thus presenting characteristics of a deeply embedded object. If this corresponds to a genuine linear motion it would be possible to detect its expected X-ray position at the time of the COUP observations back in 2003 by extrapolation, however, in the X-ray COUP images there is no detection at the extrapolated position in 1" around apart from the source COUP 510 and the next closest X-ray source is COUP 533 at 3. 1 away. Therefore a linear motions is not a plausible scenario for this source's motion. We then looked into the VLA archive for high angular resolution multiepoch observations towards this source and found 6 additional observations. These observations have been previously reported by Gómez et al. (2005Gómez et al. ( , 2008, and we followed their data calibration and imaging procedures. The phase centers of these observations are consistent with the one used in this work, except for two observations with offsets of 30 (1991.68) and 45 (2004.84), which are still smaller than the primary beam size at the observed frequency. The data calibration follows the standard procedure recommended for the pre-upgrade VLA, and the imaging discarded visibilities provided by baselines longer than 100 k to diminish the noise caused by the poorly mapped extended emission (for details see Gómez et al. 2008). The positions were corrected to consider the updated position of the phase calibrator. In total, the    Table 4. radio source is detected in 8 epochs above 5 , with rms noise levels ranging from 3 − 80 Jy bm −1 , spanning almost 25 years including the epochs 2012 and 2016 used in this work. The main image parameters of these observations are listed in Table 5 together with the measured positions of the source. Figure 10 shows the PM diagram in RA (left) and DEC (right) of the source 117 using the multiepoch observations listed in Table  5. We include the additional VLBA observation from 2015 Dzib et al. 2020). These additional data show that the linear PM estimated only from the two VLA epoch 2012 and 2016 is not representative and may indeed correspond to different components of what seems to be a binary or higher order stellar system. Proper motion measurements for this source have been previously reported in Dzib et al. (2017) using multi-epoch VLA observations with a time baseline of ∼29 years. The reported PM in their work are ( ) = 1.9 ± 1.5 mas yr −1 and = 4.2 ± 5.6 mas yr −1 which is compatible with a non moving source. The VLBA detection included in this work is compatible with the detection from VLA 2016. What is still intriguing is the fact that in every single observation only one component was detected and never both or more of them at a time, and in four visits with VLBA it was only detected once. Orbital details of this system thus remain unclear.

Radio Variability
The study presented in Forbrich et al. (2017) using the deep VLA observations showed evidence of variability at very short timescales in the range of minutes with a few cases showing changes in flux density by a factor of 10 in less than 30 minutes. Based on these results we re-imaged the central pointing of our data into time slice images of ∼ 5 minutes integration time following the procedure described in section §2 and then produced light curves for all the 272 sources in the nominal catalog of the central pointing by extracting their flux information from each of the individual 5-minute images using the methodology described in section §3.1. This time resolution also include exactly 3 science scans to ensure an even time on target throughout the time series. It resulted in a total of 41 individual images with a mean rms noise of 29 Jy bm −1 ranging between 26 − 36 Jy bm −1 . The first important issue here is the increase in the mean rms noise compared with the averaged ∼ 5 h image on where the overall rms noise is 10.24 Jy bm −1 (see Table  1). Regardless of the decrease in sensitivity all the sources were detected in at least one of the individual images, 46% of them were detected in at least half of the 41 images and 22% of them were detected in every single image.
Prior to the quantification of flux variation it was necessary to look at any systematic fluctuation. A very good test-case is the high-mass young stellar object BN, a non-variable thermal radio emitter (Forbrich et al. 2008(Forbrich et al. , 2016 with a peak flux density in the averaged 5 hrs image of = 2.370 ± 0.023 mJy bm −1 . Its peak flux density from the light curves has a mean value of = 2.354 ± 0.007 mJy bm −1 with a relative standard deviation of ∼ 5% (0.127 mJy bm −1 ). This systematic variation is smaller than the typical value for even just moderate variable sources with relative standard deviation >20% and therefore it does not impact our results of variability analysis. This measure for the BN source is also consistent with its peak flux density observed in the deep VLA catalog of = 2.3413 ± 0.0026 mJy bm −1 which represent a variation of just 0.9 ± 0.7%.
A few examples of light curves are shown in Figure 11 for sources presenting a small flare-like event and others presenting either increasing or decreasing light curves. The peak flux density of a source in the averaged long exposure image evidently does not always represent the overall evolution of its radio emission and there is a clear fluctuation at short time scales of just minutes that even not being extreme events can still be considerably brighter than the averaged peak flux density. Following Forbrich et al. (2017), an extreme variability event refers to a change in flux density by a factor of at least 10 in less than an hour. We define a variability factor VF as the ratio between the maximum and the minimum peak flux density in the timeseries dataset. Although, there is evidence of variability in several sources in the sample, none of them show changes in flux density above a factor of 5 and indeed ∼85% have VF<3.
The occurrence rate of extreme radio variability in the ONC has been estimated by Forbrich et al. (2017)     observations where they found a mean time between these events of 2220±1280 h (∼ 3 months). The sources involved in this estimate are only sources with X-ray counterpart in the COUP survey assuming that this is the most complete sample of YSOs in the inner ONC. In our sample only 157 out of 272 sources have X-ray counterparts in the COUP survey, and the cumulative radio observing time of these sources is 785 h (5 h on each of the 157 sources) for which we found no extreme variability (not by a factor 10). However, this result lies within the probabilities of only one such event in a cumulative radio observing time of 2220±1280 h and yet, our cumulative radio observing time is still below the lower limit of 940 h considering the error. The new cumulative radio observing time adds up to 7445 h including the deep (Forbrich et al. 2016) and the new observation (only the central pointing), it leads to a mean time between extreme radio variability events of 2482±1433 h. This larger uncertainty compared to the previous estimate is a result of the same number of events we are still considering (3 extreme events from Forbrich et al. 2017).
On the other hand, it is also possible to look for variability on longer timescales by comparing our central pointing against the deep VLA observations. This is a particularly good comparison since these are identical observations and all the sources are affected by the same primary beam response correction due to their equivalent distance to the phase center. In our sample 253 sources are detected in the two datasets. We can measure the variability factor between these two epochs for each source by taking the ratio between their nominal peak flux density from both catalogs. There are two sources presenting a VF>10 at this long timescale and a third one marginally below this limit with a VF of 9.4±0.6. These three sources have Xray counterparts in COUP with hydrogen column densities in the range log( ) = 21.07 − 21.52 cm −2 (Getman et al. 2005b). These parameters are summarized in Table 6.
Non-detections in the new observations of previously detected sources are not just due to a difference in sensitivity but potentially also due to variability at least for the brightest sources in the deep observations that should be clearly detected in the new VLA data. There are in total 303 sources in the deep catalog not detected in the central pointing of the new VLA data. We looked for their positions in the new data and considered 5 times the local rms within a box of 14 pixels (1. 4) as upper limits to compare against the peak flux density in the deep catalog. This allows us to estimate a lower limit for their variability factor. As expected, most of the sources have a VF below 10 and are likely too faint for the new observations, but the are two interesting sources with VF of 18 and 29 listed in Table 7. These two sources are indeed in the list of only 13 extremely variable sources in Forbrich et al. (2017). Source [FRM2016] 422 in Forbrich et al. (2017) shows an extreme radio flare on timescales of ∼40 h while source [FRM2016] 515 shows the most extreme variability within just 30 minutes (by a factor of ∼100), also coinciding with an almost simultaneous X-ray flare of similar duration. While source [FRM2016] 422 was also not detected in any of the outer pointings, source [FRM2016] 515 was detected in pointings 3 and 6 with peak flux densities of 0.15±0.01 mJy bm −1 and 0.18±0.01 mJy bm −1 , respectively, leading to a long term variability with VF∼7 respect to Forbrich et al. (2017), and VF∼5 respect to the central pointing (considering the 5 upper limit) which were observed with one month difference. Both sources [FRM2016] 422 and [FRM2016] 515 have reported spectral types K8 and O9.5-B2, respectively (Hillenbrand et al. 2013), although for the high-mass star it is noted that the radio emission may be detected from unresolved lower-mass companion.
Finally, another interesting finding concerns the new detections in the central pointing not detected in the previous deep VLA data. There are 19 new detections, 10 of which have X-ray counterparts in the COUP survey. There are 9 sources with S/N>10 and two of them have high S/N of 50 and 82. The main parameters of these two bright new sources are listed in Table 8.

SUMMARY AND CONCLUSIONS
We have presented a new deep, high-resolution catalog of compact radio sources towards the ONC at centimeter wavelengths using the most extended configuration of the VLA. This is the deepest catalog for the surrounding areas of the ONC reported to date, reaching rms noise levels between 3−5 Jy bm −1 at distances of ∼10. 4 from the center of the cluster, significantly improving the general census of known compact radio sources in the cluster. We detected a total of 521 radio sources above 5 threshold over an area of ∼20 × 20 . In this catalog, 198 sources are new detections not previously reported at these frequencies. The highest stellar surface density occurs in the inner region, and yet the number of radio sources may still be underestimated due to the difficulty of disentangling small-scale structure of the nebula and stellar point sources.
With our new catalog, we are sensitive not only to stellar radio emission but also to radio emission originating elsewhere in the ONC, for example in outflows and shocks. It turned out in this regard that even the relatively short time baseline of 4.19 years is sufficient to trace high-velocity proper motions in the ONC, which can be more easily measured with phase-referenced radio than with optical observations. We identify radio sources in the deep 2012 catalog as co-moving with and thus originating in ejecta of the OMC1 explosion. While we also appear to detect fast proper motions towards stellar sources, those are likely due to insufficiently sampled multiple systems.
Column (2): VF between the two epochs considering 5 times the local rms in the new observations. Column (5) and (6)   estingly, the surface density of radio sources falls off faster than that of X-ray sources, tracing young stars. This may be due to a surplus of thermal/free-free sources in the inner ONC that are ionized by 1 Ori C. Additionally, we find that most (> 50%) radio sources have X-ray counterparts, throughout the area studied here, while only a minority (< 20%) of X-ray sources have radio counterparts. X-ray emission thus remains a poor predictor of radio emission in this sample of young stars. Given the regional differences that we find, this is not simply a sensitivity effect. Finally, a radio variability analysis for sources in the inner ONC is presented as a follow-up of our previous deep VLA variability study. We produced radio light curves at high time resolution for sources in the central pointing finding changes in flux density by a factor 5 on timscales of minutes to a few hours thus we do not find extreme variability events. The majority of the sources have detections in the previous deep VLA observations enabling the study of long term variability where we only find two sources with changes in flux density 10. Based in these two studies we find a mean time between extreme radio variability events of 2482± 1433 h.