Very Long Baseline Interferometry Observations of the Proposed Radio Counterpart of an EGRET Source

We present high-resolution radio interferometric imaging observations of the radio source NVSS J182659+343113 (hereafter J1826+3431), the proposed radio counterpart of the ${\gamma}$-ray source, 3EG J1824+3441 detected by the Energetic Gamma Ray Experiment Telescope (EGRET) on board the Compton Gamma Ray Observatory satellite. We analyzed eight epochs of archival multi-frequency very long baseline interferometry data. We imaged the asymmetric core-jet structure of the source, and detected apparent superluminal motion in the jet. At the highest observing frequency, 15.3 GHz, the core shows high brightness temperature indicating Doppler boosting. Additionally, the radio features undergo substantial flux density variability. These findings strengthen the previous claim of the association of the blazar J1826+3431 with the possible ${\gamma}$-ray source, 3EG J1824+3441.


Introduction
Blazars are radio-loud active galactic nuclei (AGN) seen at a very small angle to the jet direction. The relativistic speed of the emitting plasma in the jet and the small angle to the line of sight give rise to various relativistic effects such as apparent superluminal motion and relativistic Doppler boosting of the radiation. Due to the latter, the emission of the approaching jet is enhanced while that of the receding jet (also called counterjet) is diminished significantly. The jet-to-counterjet flux density ratio can be as high as 1000 [1]. As a consequence, only the approaching one of the originally symmetrically launched jets can be observed, thus blazars appear asymmetric in that regard. The jet emission is explained as synchrotron radiation from electrons and/or positrons moving with relativistic speeds in magnetic fields. The jet synchrotron radiation is responsible for the lower-energy bump seen in the broad-band spectral energy distribution (SED) of blazars, peaking between millimetre-wavelength and X-ray regimes [2].
Blazars constitute the most populous group of extragalactic γ-ray emitter sources according to the observations of the Fermi Large Area Telescope (LAT; e.g., [3]). The γ-ray emission, responsible for the higher-energy bump in the blazar SEDs, is usually attributed to inverse-Compton process in leptonic models (e.g., [2] and references therein) or to synchrotron radiation from heavy charged particles, e.g., protons ( [4] and references therein) in hadronic models. The seed photons of the inverse-Compton

VLBI Data
In the very long baseline interferometry (VLBI) image database of the astrogeo.org website (http://astrogeo.org/ maintained by L. Petrov), there are four epochs of dual-band 2.3-and 8.3/8.7-GHz observations and three epochs of 15.3 GHz measurements available of J1826+3431 between 1996 and 2018. In each of these observations, all ten antennas of the VLBA participated. Further details of the observations are given in Table 1. Additionally, the source was observed with the European VLBI Network (EVN) at 1.7 GHz on 21 February 2014 (project code: EF025, PI: S. Frey). This observation was conducted in phase-referencing mode [14], and J1826+3431 was used as the phase-reference calibrator source for one of the faint targets in that project [15]. Eight antennas of the EVN participated in the measurements: Effelsberg in Germany, the Jodrell Bank Lovell Telescope in the United Kingdom, Medicina and Noto in Italy, Onsala in Sweden, Toruń in Poland, the Westerbork Synthesis Radio Telescope in the Netherlands, and Sheshan in China. Further details of the observation are given in Table 1 and in [15].

VLA Data
The public data archive (https://archive.nrao.edu/) of the U.S. National Radio Astronomy Observatory (NRAO) lists various data sets which contain VLA observations of J1826+3431. To obtain supplementary information on the arcsec-scale extended radio structure, as well as the possible variability of the source, we selected two data sets obtained at 4.86 GHz for analysis. The experiment AS291 (PI: W.C. Saslaw) was observed at the epoch 1987.493 in the most extended A configuration of the array that provides the highest angular resolution, in our case, ∼0.5 arcsec. The second experiment from 2009.388 (project code: CALSUR) used the less extended VLA B configuration, providing resolution of a few arcsecs. The integration times spent on J1826+3431 were 219 s and 43 s in 1987 and 2009, respectively. Both experiments used 100 MHz bandwidth.

VLBA Data
Calibrated visibility files were downloaded from the astrogeo.org website, and their hybrid mapping procedure was performed using the Caltech DIFMAP [16] software package. Asterisks (*) denote those epochs in Tables A1-A4 where we time-averaged the calibrated data into 10-s blocks. In a process of hybrid mapping, the CLEAN algorithm [17] was used and phase-only self-calibration was done. At the end of the procedure, steps of amplitude and phase self-calibration were performed, gradually for shorter and shorter solution time intervals down to a few minutes. Images ( Figure 1) were created using natural weighting, with the visibility errors raised to the power −1.

EVN Data
The EVN data calibration was done in the usual manner (e.g., [18]) using the NRAO Astronomical Image Processing System (AIPS, [19]). The detailed description is given in [15]. The calibrated data set was imaged in DIFMAP [16] using the hybrid mapping technique, which involved several steps of CLEAN iterations interleaved by phase-only self-calibrations. When the signal-to-noise ratio in the residual image could not be improved any further, amplitude and phase self-calibration were performed. First, it was done for the whole observation; afterwards, the corrections were computed for a subset of antennas and with gradually decreasing solution time intervals down to 10 min. The image obtained is displayed in Figure 1a.

VLA Data
We performed standard amplitude and phase calibration in AIPS. The VLA flux density calibrator to set the amplitude scale was 3C 48 in both experiments. The calibrated visibility data were exported from AIPS to DIFMAP for hybrid mapping. The resulting CLEAN images are shown in Figure 2.

Model Fitting
To describe the brightness distribution for a quantitative analysis of the source structure, we applied the technique of model fitting directly to the calibrated interferometric visibility data [20], for both the VLBI and VLA measurements. We used DIFMAP to represent the source brightness distribution with a set of circular Gaussian model components.      Table 2. Details of the contour maps shown in Figure 1. The restoring beam major axis position angle (PA) is measured from north through east. Further positive contour levels in Figure 1 increase by a factor of 2.

Results
Based on VLBI data, the inner structure of J1826+3431 shows a southeast oriented core-jet morphology consistently at all frequencies and epochs (see Figure 1 for selected example maps). The core and typically 4-5 jet features were needed to adequately model the source with circular Gaussian components at 2.3, 8.3/8.7, and 15.3 GHz. At the lowest frequency, 1.7 GHz, only two model components were necessary to describe the jet brightness distribution. During the fitting procedure, all of the parameters, i.e., positions, flux densities, and full width at half-maximum (FWHM) sizes of the components were allowed to vary, except for one epoch. At 2005.427, the fitting procedure was not able to constrain the FWHM size of the core component; therefore it was fixed at the smallest angular size resolvable with the interferometer, 0.04 mas. To calculate this value, we used the formula given by [21]. Thus, we can only give an upper limit of the size of the core component at that epoch. The parameters of the best-fit models are summarized in Tables A1-A4; the model component positions are also denoted in the images in Figure 1.
When deriving the uncertainties of the fitted parameters of the VLBI observations, we assumed 10% flux density calibration uncertainties following [22]. The positional uncertainties were calculated conservatively following [23]. We quadratically added the uncertainties determined using the relations given by [24] to the 10% of the restoring beam size in the cases of the core components, and to the 20% of the restoring beam size when deriving the positional uncertainty of less bright jet components. To calculate the uncertainties of the FWHM sizes, we used the formula of [24].
The arcsec-scale radio structure of J1826+3431 revealed by our 4.86-GHz VLA images ( Figure 2) indicates a continuing jet towards the southeast, in approximately the same position angle as the mas-scale VLBI jet ( Figure 1). The jet can be traced out to ∼2.5 (corresponding to ∼20, kpc projected distance from the core) in the higher-resolution A-configuration image. There is an indication of a sharp right-angle bend towards south, as supported by the lower-resolution B-configuration VLA image where the radio emission seems extended in that direction ( Figure 2). The sum of the flux densities in the fitted Gaussian model components (428 ± 21 mJy at 1987.493 and 333 ± 17 mJy at 2009.388) provide evidence for significant variability.

Component Proper Motions in the Mas-Scale Jet of J1826+3431
To identify the jet components throughout the full time span of the VLBI observations, we considered their positions, FWHM sizes, and flux densities. The components that are identified across the different epochs are denoted by the same identifiers in Tables A2-A4. The core component is numbered with 0 at each epoch and frequency; thus, it is denoted as L0, C0, D0, and J0 at 1.7, 2.3, 8.3/8.7, and 15.3 GHz, respectively.
For further analysis, we assumed the core to be stationary and we derived the positions of the fitted jet components relative to it (Tables A1-A4). The jet components are assumed to move away from the core at constant speeds along a straight trajectory. Thus, we fitted their core separations with linear functions. The derived angular proper motions were converted to apparent speeds in units of the speed of light, c, using the relation [1] where µ is the angular proper motion of the component in rad s −1 , d L is the luminosity distance in m, c is the speed of light in m s −1 , and z is the redshift of the source. At 2.3 GHz, at most five components were needed to describe the jet structure. The inner ≤10 mas region of the jet could not be resolved into two distinct components in 1996.370, unlike at later epochs (components C1 and C2). Therefore, the single fitted inner jet component, Ca, cannot be unambiguously identified with either of the components detected later. Thus, the motion of C1 and C2 can only be tracked at the relatively closely-spaced observing epochs of the 2010s. Because of the short time baseline, these two components show no significant change in their core separation during the available observing epochs. On the other hand, component C3 could be unambiguously identified through all epochs. It shows an angular proper motion of µ = 0.12 ± 0.07 mas yr −1 , which corresponds to an apparent speed of β app = 9.4 ± 5.5. At large core separations, component C4 could only be detected at two epochs. Therefore, it is unsuitable for further proper motion calculations. The core separations for each component and the fitted linear proper motion curve for C3 are shown in Figure 3. At 8.3/8.7 GHz, three jet components, D2, D3, and D4, can be detected at all four epochs, while D1 can only be detected at the last three epochs. This may indicate that it was ejected some time between the epochs 1996.370 and 2015.064. No discernible proper motions could be detected for components D1 and D2, while D3 and D4 show apparent superluminal motions ( Table 3). The core separations and the fitted linear curves for D3 and D4 are shown in Figure 3. As we performed proper motion and jet parameter calculations using D3 and D4, we present contour maps of J1826+3431 made at this frequency at all four epochs in Figure 4, indicating the fitted model components.   Table 4.  Table 4. Details of the contour maps shown in Figure 4. The restoring beam major axis position angle (PA) is measured from north through east. Further positive contour levels in Figure 4 increase by a factor of 2. To summarize, apparent superluminal motions could be detected in three components at relatively larger core separations (C3, D3 and D4) in the jet of J1826+3431. Interestingly, the highest apparent speed is shown by a component (D4) located at the largest core separation (∼10 mas) implying acceleration in the jet. Acceleration is often detected in the jets of blazars (e.g., [25]), and although it usually takes place at smaller core separations, there are examples for faster components further down the jets (e.g., PMN J0405−1308 [25,26], PKS 2201+171 [27]).

Relativistic Beaming in J1826+3431
At the highest observing frequency, 15.3 GHz, which provides the best angular resolution, we derived the brightness temperature of the core component using the following equation [28]: where S is the flux density of the fitted component in Jy, θ is the FWHM diameter of the circular Gaussian in mas, and ν is the observing frequency in GHz. In 2003.655, we were unable to resolve the core region into two components as at the later epochs. As a test, we tried to fit the innermost region with an elliptical Gaussian component and the resulting fit confirmed that the core was elongated in the direction of component J1 detected at the later epochs. Therefore, J0 and J1 were blended in 2003.655, thus the FWHM size and the flux density are overestimates of the core parameters at this epoch. Since the brightness temperature is inversely proportional to the square of the FWHM size, its derived value can be regarded as a lower limit, ≥ 2.8 × 10 11 K. In 2004.616, the obtained brightness temperature is T b = (12.3 ± 1.4) × 10 11 K. At epoch 2005.427, we could only derive an upper limit for the size of the core component, thus only a lower limit for the brightness temperature can be given, ≥12 × 10 11 K. The lower limits agree with the T b estimate made at the middle epoch within its errors.
All brightness temperature values exceed the equipartition limit, T int ≈ 5 · 10 10 K [29], indicating relativistic beaming in the source. The relativistic beaming can be quantified by the Doppler boosting factor, δ vlbi , which can be estimated from the measured T b value and the equipartition brightness temperature limit using the following relation [1]: The obtained brightness temperature value corresponds to a Doppler factor of δ vlbi = 24.6 ± 2.9.

Jet Parameters
The jet parameters, the bulk Lorentz factor (γ), and the inclination angle with respect to the line of sight (φ) can be calculated from the Doppler factor and the apparent jet speed as [1] and tan φ = 2β app The component close to the innermost core region with the best-defined positions and thus proper motion is D3; therefore, we used the apparent superluminal speed obtained for this jet feature, β D3 app = 4.8 ± 1.6, to derive the jet parameters. Using the Doppler factor obtained from the 15.3-GHz measurement at epoch 2004.616, we get a Lorentz factor of γ = 12.8 ± 1.4, implying relativistic bulk jet velocity in units of the speed of light, β = 0.997 ± 0.002. The corresponding jet inclination angle is φ = 0.9 • ± 0.3 • .

Radio Spectrum of J1826+3431
We collected the radio flux density measurements of J1826+3431 from the literature ( Table 5). We fitted a power-law curve to these data points to obtain the spectral index (α; defined as S ∝ ν α , where S is the flux density and ν is the observing frequency). Since the flux density measurements are not simultaneous, source variability may affect the spectral index calculation, so the value can be considered as tentative. The resulting fit is shown with a black line in Figure 5, the obtained spectral index is α = −0.19 ± 0.05, indicating a flat spectrum.
For comparison, we also show the sum of the VLBI Gaussian model component flux densities with blue squares and upward arrows in Figure 5. For each observing frequency, we have chosen the observation that is the closest to the 1.7-GHz EVN measurement. VLBI observations are insensitive to the large (arcsec) scale radio structure; therefore, in the presence of such emission, they underestimate the total flux density of a source. This can be seen in Figure 5 where all VLBI points are indeed below the fitted spectrum. Archival data VLA-B VLA-A VLBI Figure 5. Radio spectrum of J1826+3431. Black circles are from archival measurements, their details are given in Table 5. Green triangles indicate the sum of modelfit components of the VLA measurements presented in this paper. Black line represents the power-law fit to the above total flux density points. Blue squares with arrows show the sum of modelfit components of high-resolution VLBI measurements from this paper, as lower limits to the total flux density. The simultaneous 2.3-and 8.3/8.7-GHz observations allow us to map the spectral index distribution of the mas-scale radio emission in J1826+3431 using VIMAP [36]. We present a spectral index image in Figure 6. We used our data from the epoch 2015.064. After using elliptical masks to exclude the cores from the images, we matched the 2.3-and 8.7-GHz maps of J1826+3431 by using two-dimensional correlation. Considering that the frequencies are relatively close to each other and that the measurements were performed at the same time, it was not necessary to shift the brightness peak. The spectral index image ( Figure 6) suggests a flat-spectrum radio emission (α ≈ −0.1 ± 0.1) in the core region and along the jet closer to the core. It is in good agreement with the spectral index calculated using the archival total flux density data ( Figure 5) since the dominant emission feature of the source is the core. The spectrum gradually steepens further away from the core along the jet, as expected from optically thin synchrotron radio emission.

Flux Density Variability
To characterize variability in the VLBI-detected features, we plot the sum of flux densities in the fitted Gaussian model components at different observing frequencies as a function of time in Figure 7. As can be clearly seen, especially for the 2.

γ-Ray Properties of 3EG J1824+3441
3EG J1824+3441 was detected with a weak γ-ray flux of (28.7 ± 9.3) × 10 −8 photon s −2 cm −1 by EGRET [6]. However, an alternate version of the EGRET catalogue [37] (see also [7]) which used reprocessed data assuming a different model for the Galactic interstellar emission does not contain as many as 107 sources previously detected in [6], including also 3EG J1824+3441. In that sense, it is not surprising that no γ-ray source is listed in the most recent fourth Fermi/LAT catalogue [3] within a 1 • radius of the position of the putative EGRET γ-ray source.
Alternatively, even if the source exists, long-term variability of its γ-ray flux can also explain the Fermi non-detection. Indeed, a comparison of high-confidence (>4σ) EGRET detections with the third Fermi/LAT catalogue revealed 10 extragalactic γ-ray sources missing from the latter [38]. Subsequent analysis of seven years of Fermi data showed that five sources were in a low-luminosity state during the Fermi observations compared to their luminosity during the EGRET observations [38]. The generally decreasing long-term trend we see in the radio flux density of the possible counterpart, the blazar J1826+3431, with respect to the 1990s when CGRO was operational (see e.g., our VLA results and Figure 7) is also consistent with a γ-ray flux decrease.

Conclusions
We analyzed eight epochs of high-resolution multi-frequency VLBI data of the radio source J1826+3431, which is the proposed radio counterpart of the EGRET γ-ray source, 3EG J1824+3441. The archival observations span a period of more than 22 years. Each data set was imaged, and the brightness distribution of the source was modeled using circular Gaussian components, in order to characterize the physical and geometric properties of the jet and to analyze the kinematics of J1826+3431. One and two moving jet components could be securely identified at multiple epochs in the 2.3and 8.3/8.7-GHz data, respectively. All of them showed outward motion with apparent superluminal speed. The estimated brightness temperature of the most compact component, the core at 15.3 GHz suggests that the radio emission is relativistically beamed in the source. From the relativistic effects seen in the radio jet and the strong variability of the core component at the measured frequencies, the radio source J1826+3431 can be identified as a blazar. Since the vast majority of extragalactic γ-ray sources are associated with blazars, we can conclude that J1826+3431 is very likely the radio counterpart of the unidentified EGRET source 3EG J1824+3441. However, no γ-ray source is found at this position in the latest Fermi/LAT catalogue [3], and its detection can be regarded uncertain based on an alternative version of the EGRET catalogue [37]. Acknowledgments: The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. The European VLBI Network is a joint facility of independent European, African, Asian, and North American radio astronomy institutes. Scientific results from data presented in this publication are derived from the following EVN project code: EF025. We acknowledge the use of archival calibrated VLBI data from the Astrogeo Center database maintained by Leonid Petrov.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Abbreviations
The following abbreviations are used in this manuscript:

Appendix A
Here, we provide the results of the model fitting to the VLBI visibility data. In Tables A1-A4, we list the best-fit parameters of the circular Gaussian components used to describe the brightness distribution of the 1.7-, 2.3-, 8.3/8.7-, and 15.3-GHz calibrated data of J1826+3431, respectively.  C4  Table A3. Parameters of the fitted Gaussian model components at 8.3/8.7 GHz. The observational epochs are given in column 1. The flux densities, relative positions in right ascension and declination directions with respect to the core component, and the FWHM sizes are listed in columns 2, 3, 4, and 5, respectively. In the final column, we give the component identifiers. Asterisks (*) denote those epochs where we time-averaged the calibrated data into 10-s blocks.

Epoch
Flux