Precessing Jet in the High-Redshift Blazar J0017+8135

The prominent flat-spectrum radio quasar J0017+8135 (S5 0014+81) at z = 3.366 is one of the most luminous active galactic nuclei (AGN) known. Its milliarcsecond-scale radio jet structure has been studied with very long baseline interferometry (VLBI) since the 1980s. The quasar was selected as one of the original defining objects of the International Celestial Reference Frame, but left out from its current second realization (ICRF2) because of systematic long-term positional variations. Here we analyse archival 8.6- and 2.3-GHz VLBI imaging data collected at nearly 100 different epochs during more than 20 years, to obtain information about the kinematics of jet components. Because of the cosmological time dilation, extensive VLBI monitoring data are essential to reveal changes in the jet structure of high-redshift AGN. In the case of J0017+8135, the data can be described with a simple kinematic model of jet precession with a 12-year periodicity in the observer's frame.


Introduction
The blazar J0017+8135 (S5 0014+81) is listed in the S5 sample of strong northern radio sources with S = 551 mJy total flux density at 5 GHz [1]. Its spectroscopic redshift is z = 3.366 [2]. This flat-spectrum radio quasar is known as one of the most luminous active galactic nuclei (AGN) and has been subject to various studies at multiple wavebands, from the radio to X-rays (e.g., [3][4][5][6]). The mass of its central black hole is M = 4 × 10 10 M , among the highest values measured in AGN [5]. The source has been extensively studied with very long baseline interferometry (VLBI) since the 1980s, and it was also one of the sources first observed with space VLBI [7,8]. The one-sided core-jet structure of J0017+8135 extending towards the south at milliarcsecond (mas) angular scales is typical for blazars. The jet components show no or moderate apparent superluminal motion [9,10]. J0017+8135 was selected as one of the original defining objects of the International Celestial Reference Frame (ICRF), but left out from its current second realization (ICRF2) because of systematic long-term positional variations [11]. In this paper, we make use of the publicly available archival 2.3-and 8.6-GHz VLBI imaging data spanning more than two decades, to investigate the positional variations of the jet components in J0017+8135. We find indication of periodic changes and attempt to explain them in the framework of a simple jet precession model. A detailed jet kinematic study for a blazar at such a high redshift is unique. Changes in the source structure appear slower by a factor of (1 + z) in the observer's frame due to cosmological time dilation; therefore, long-term VLBI monitoring is necessary to reliably reveal them in the case of high-redshift radio-loud AGN.

VLBI Observations, Images and Brightness Distribution Modeling
The study presented here is based on 95 VLBI imaging data sets at 8.6 GHz (X band). In 74 cases, VLBI data obtained at the same epochs are also available at 2.3 GHz (S band). The observations were done with the U.S. National Radio Astronomy Observatory (NRAO) Very Long Baseline Array (VLBA), occasionally supplemented with various other radio telescopes, covering the time range from August 1994 to February 2015. We used the calibrated visibility data sets downloaded from [12]. The imaging and model-fitting were performed through standard procedures using Difmap [13] in an automated way. The visibility data were modeled with sets of 2 or 3 Gaussian brightness distributions. The core was fitted with an elliptical Gaussian component at both frequencies. The jet at 8.6 GHz was fitted with one, and at 2.3 GHz with two circular Gaussians. A pair of typical VLBI images of the source with the fitted model components is presented in Figure 1. The core is labeled as C, the jet components are J1, J2 and J3, respectively. In the 2.3-GHz images, the central component is labeled as C+J1 as it blends the two innermost components found in the higher-resolution 8.6-GHz data. Although an opacity correction would be needed, it would have a minor effect due to the relatively close observing frequencies. Thus, in both cases, the central component was considered as the kinematic centre for the modeling.

A Simple Jet Precession Model
We fitted a kinematic model [14] to the relative positions of the brightness distribution components measured as a function of time, to quantitatively describe the periodic behaviour of their motion. We assumed that the plasma blobs ejected at different epochs move in a ballistic trajectory. At each epoch, we see the ejected components at their actual projected positions, and the base of the jet as the optically thick 'core' region. The relative positions of the identified components can change over time if the core is wobbling (i.e., changing its absolute position). Thus, if we choose the core component as the kinematic centre of the jet, the movements of the core will be reflected in the relative movements of the ejected components. The relative positions of the jet components show periodic difference from the ballistic track, thus we assume a possible precession of the jet. In our precession model, we make three assumptions for the kinematics and geometry of the jet: (1) the jet sweeps a conical surface whose half opening angle Ω and axis are constant over time; (2) the jet is precessing with a constant ω angular velocity around the cone axis; and (3) the jet components ejected at t 0 epoch can be described at t with two time-dependent angles: η(t) and Φ(t). In this scenario, the ejected components will follow a helical path in the plane of the sky. The projections of this kind of motion to the right ascension (RA) and declination (Dec) axes can be described as superpositions of a linear and a periodic function. The angular frequency of the projected motions is characteristic to the jet precession, thus it has to be equal in both directions. The apparent precessing motion of the jet is the combination of the two perpendicular oscillations in this model. Figure 2 shows the geometric configuration of the jet model. Assuming this configuration, for each jet feature, we can fit a superposition of a linear and a periodic function to the projected component RA and Dec angular distances from the core: where i ∈ [RA, Dec], s i (t) is the component distance from the core at the time t, s 0i is the component distance from the core at t = 0, µ i is the proper motion of the component, ω is the common angular frequency, A i and ϕ i are the fitted amplitude and phase, respectively. As a boundary condition, the model requires that the relative RA and Dec positions of the components at the ejection epoch t 0 are However, if the amplitude of the precession is small, we can simplify this condition so that only the fitted ballistic slopes have to intersect the t axis at the same point. From this assumption, we get Thus we have to fit the model to the relative RA and Dec position data for all components at the same time, with the common parameter ω. Considering the simplified boundary condition, the kinematic model can be expressed as

Jet Kinematic Model Parameter Estimates
To analyze the precession of the jet model components, we fitted Equation (4) to the jet component positions relative to the core as derived from the multi-epoch VLBI data. We applied the Metropolis-Hastings algorithm [15,16] to look for the best-fit parameters for each component.  Based on the best-fitting model parameters, we can express the main kinematic properties of the jet components (Table 1). We estimated the angular frequency ω, the precession period T, and the apparent proper motions β app for each component (expressed in the units of the speed of light c). To estimate φ(t 0 ) and Ω with reasonably small error, we would have to observe more model components in the jet [14]. Nevertheless, we also estimated the projected angle of the precession cone axis on the plane of the sky at the ejection epoch, η(t 0 ), measured from north through east. ω (year −1 ) 0.51 ± 0.08 0.51 ± 0.15 0.5 ± 0.1 T (year) 12 ± 2 12 ± 4 12 ± 2 β app 0.9 ± 0.2 180 ± 15 195 ± 34 189 ± 10

Discussion
We found that, according to the expectations, the ejection epochs were in chronological order for the J3, J2, and J1 components, respectively. The apparent motion of the components is dominantly in the Dec direction, and the amplitude of the periodic motion is higher in this direction. Thus the precession cone is projected onto the sky as an ellipse stretched in the north-south direction, in accordance with the absolute astrometric position changes [11]. Somewhat similarly to earlier results based on VLBI observations spanning shorter intervals [9,10], our analysis reveals mild superluminal motion (∼1 − 2c) in the jet components of J0017+8135 but with smaller uncertainties.
In the framework of the kinematic model applied here, the jet has a precession period of (12 ± 3) year in the observer's frame. This is a unique result for a blazar at z > 3 where the cosmological time dilation makes it more difficult to detect periodicities compared to low-redshift sources. However, it should be noted that this analysis is based on data that span an interval shorter than twice the putative period. Therefore any stochastic variability might be mistakenly identified as periodic. Continuing VLBI monitoring of J0017+8135 on a longer term, and the analysis of absolute astrometric position variations [11] and total flux density monitoring data [17] in the context of the jet precession model would help confirming our results in the future.
Finally, it is also possible that physical mechanisms other than precession cause the observed apparently helical path of the jet components. For example, plasma instabilities propagating along a conical jet can cause this effect. In fact, studies of large samples of sources monitored with VLBI revealed that non-ballistic motion is ubiquitous in pc-scale AGN jets ( [18] and references therein).