A parsec-scale wobbling jet in the high-synchrotron peaked blazar PG 1553+113

PG 1553+113 is the first blazar showing an approximately two-year quasi-periodic pattern in its gamma-ray light curve. Such quasi-periodicity might have a geometrical origin, possibly related to the precessing nature of the jet, or could be intrinsic to the source and related to pulsational accretion flow instabilities. By means of a ~2yr very long baseline array (VLBA) monitoring at 15, 24, and 43 GHz we investigate the source pc-scale properties during an entire cycle of gamma-ray activity in the period 2015-2017. In contrast to the well-defined periodicity in the gamma-ray emission, at radio frequencies no clear periodic pattern can be recognized. The jet position angle, constrained by means of the total intensity ridge line, varies across the different observing epochs in the range 40-60 deg. We also investigate the time evolution of the source polarization properties, including the rotation measure. The brightness temperature is found to decrease as the frequency increases with an intrinsic value of ~1.5 x 10^10 K and the estimated Doppler factor is ~1.4.


Introduction
Blazars are the most extreme objects in the family of active galactic nuclei (AGNs). With their jets closely aligned to our line of sight, they represent the best target for the study of both the physics of particle acceleration and the role played by magnetic fields in these extreme plasma environments (e.g. Blandford & Königl 1979;Marscher et al. 2008). Blazars include both flat spectrum radio quasars (FSRQs -with an optical spectrum dominated by strong emission lines) and BL Lac objects (BL Lacs -showing almost featureless spectra with occasional weak optical emission lines, Stickel et al. (1991)). Blazar emission spans the entire electromagnetic spectrum from radio up to TeV γ-ray energies, and is mostly non-thermal in nature and linearly polarised, providing us with important insights into the magnetic field structure (e.g. Gabuzda et al. 1994;Gómez et al. 2008;Orienti et al. 2011;Hovatta et al. 2012;Casadio et al. 2019). Blazars tend to show erratic variability across the electromagnetic spectrum on different timescales, ranging from years (Ulrich et al. 1997) down to a few minutes (e.g. Aharonian et al. 2007; Al- Email: rlico@mpifr.de bert et al. 2007b). A small subset of blazar sources exhibit possible quasi-periodic variability across the radio, optical, and X-ray emission bands (e.g. Villata et al. 2004;Hovatta et al. 2008;Valtonen et al. 2008;Wiita 2011;King et al. 2013). Over the last few years, thanks to the Fermi Large Area Telescope (Fermi-LAT) continuous monitoring of the sky in the MeV-GeV energy range (Atwood et al. 2009), quasi-periodic variability has been investigated at γ-ray energies in a number of blazars (e.g. Prokhorov & Moraghan 2017). Such high-energy periodicity might be related to jet precession and/or modulation of the accretion rate onto the central engine(s). As such, the study of quasi-periodic variability at γ-ray energies can help shed light on fundamental issues such as the disc-jet connection and the nature of the jet's magnetic fields, and could provide further insight into gravitational wave production in binary super massive black hole (SMBH) systems (Abbott et al. 2016).
∼ 2.18 ± 0.08 years (Ackermann et al. 2015) has been observed in its γ-ray light curve, providing us with a unique opportunity to investigate this peculiar behaviour.
PG 1553+113 exhibits an almost featureless optical spectrum (Miller & Green 1983;Falomo & Treves 1990) and the emission at all wavelengths, from radio to very high energies (VHE; E> 100 GeV), can be attributed to the non-thermal jet emission, without significant contributions from other substructures (e.g. the disc, the corona, or the broad line region). The host galaxy is yet undetected and the redshift determination is still uncertain. Using a method based on Bayesian statistics and the extragalactic background light absorption effects on the VHE spectrum, Abramowski et al. (2015) obtained a value of z = 0.49 ± 0.04. Additional indirect methods, such as the Lyα forest method and the absence of a break in the VHE spectrum, constrain the redshift to the range 0.3-0.6 (e.g. Qin et al. 2018;Landoni et al. 2014;Danforth et al. 2010;Prandini et al. 2010;Mazin & Goebel 2007). This source is classified as a high synchrotron peaked (HSP) blazar (Giommi et al. 1995) and was detected by the HESS and MAGIC γ-ray Cherenkov telescopes in 2005 (Aharonian et al. 2006;Albert et al. 2007a).
The quasi-periodic trend found in the γ-ray light curve of PG 1553+113, although still debated (e.g. Covino et al. 2019), is matched by similar behaviour at lower frequencies. In addition to the periodicity in the γ-ray emission, hints of periodicities have been found in the optical emission light curve spanning ∼10 years (2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014) obtained within the Tuorla optical monitoring program 1 and in the 0.3-10 keV light curve obtained by XRT on board the Neil Gehrles Swift satellite, although at a lower significance level due to the sparse observations (Tavani et al. 2018).
To investigate the radio properties of the innermost regions of the jet, where the γ rays are thought to be produced, highresolution very long baseline interferometry (VLBI) observations are needed. On VLBI scales, the source shows a compact nuclear region with jet structure in the northeast direction (position angle ∼ 50 • ), that is well detected out to roughly 20 pc from the core region (Rector et al. 2003). Beyond this distance the jet becomes very faint and diffuse. Although several multiwavelength (MWL) observing campaigns of PG 1553+113 (e.g. Raiteri et al. 2015;Hovatta et al. 2016;Itoh et al. 2016;Raiteri et al. 2017) were conducted in recent years, only sparse VLBI observations of PG 1553+113 can be found in the literature, including the 15 GHz observations within the MOJAVE database 2 and the TeV Blazars very long baseline array (VLBA) monitoring program led by B. G. Piner and P. G. Edwards 3 . For these reasons, we conducted a systematic multi-frequency (15, 24, and 43 GHz) VLBA monitoring of PG 1553+113 covering a period between two consecutive maxima in the γ-ray activity during 2015-2017. The main purpose of this observing campaign is to investigate and characterise the parsec-scale source properties and their evolution with time as well as the possible connections with the observed flux density modulation.
In this paper we describe the observing campaign and the data analysis in Sect. 2, while the main results are presented in Sect. 3. A general discussion and summary of our findings are presented in Sects. 4 and 5, respectively. Throughout the paper we use a ΛCDM cosmology with h = 0.71, Ω m = 0.27, and Ω Λ = 0.73 (Komatsu et al. 2011). The spectral index α is defined as S ν ∝ ν −α , and all angles are measured from north to east. At a

Observations and data analysis
We monitored PG 1553+113 with the VLBA at 15, 24, and 43 GHz in full polarisation from February 2015 to June 2017. The observations are based on two related observing projects: two six-hour observing sessions during the period of its maximum γ-ray activity in 2015 (BL214) and a full two-year monitoring program with roughly bimonthly runs (BL215, 54-h in total). Table 1 reports the log of the observations. In some epochs one or more VLBA stations were missing or flagged due to technical problems (as reported in the last column in Table 1); when a station is missing only for a specific observing band, it is indicated in brackets. We were able to obtain absolute electric vector position angle (EVPA) orientations for the four observing epochs (MJD 57055,57078,57650,57755) for which we had quasisimultaneous Karl G. Jansky very large array (JVLA) observations at similar frequencies. We used the source J1310+3220 as the instrumental polarisation calibrator. We also made use of the 15 GHz observations provided by the Owens Valley Radio Observatory (OVRO) blazar monitoring program (Richards et al. 2011). PG 1553+113 has been monitored with OVRO since 2008 with a cadence of about two observations per week. In this work we use the observations over the MJD range 54696-58153 4 . We used the software package Astronomical Image Processing System (AIPS; Greisen 2003) for the data calibration, the fringe-fitting, and the detection of cross-polarised fringes. We determined the instrumental polarisation leakages (D-terms) with the AIPS task LPCAL. For the final images we used the CLEAN and self-calibration procedures in the DIFMAP software package (Shepherd 1997). The core flux densities were obtained by fitting elliptical Gaussian components to the calibrated visibilities for each epoch at each frequency in DIFMAP. In this work the core is identified as a bright and stationary feature at the upper end of the jet, with a roughly flat spectrum and high brightness temperature (of the order of 10 10 K, see Sect. 3.5). For producing the polarisation images presented in this work we made use of IDL routines developed by J. L. Gómez and for the The overlaid lowest total intensity contour is at 0.4%, 1.3%, and 1.8% of the peak at 15, 24, and 43 GHz (see Table 2), respectively, with the following contours a factor of two higher. The colour scale represents the linearly polarised intensity.
Article number, page 3 of 11 A&A proofs: manuscript no. Lico_PG1553 Table 2. Summary of the total and polarised intensity parameters shown in Fig. 2 and in the lower frame of Fig. 6 Notes.
(a) Peak flux density (S peak ) in mJy; (b) uncertainties on S peak ; (c) core region flux density (S core ) in mJy; (d) uncertainties on S core ; (e) Total flux density (S total ) in mJy; (f) uncertainties on S total ; (g) polarised flux density (P) in mJy; (h) uncertainties on P; (i) fractional polarisation (m); (j) uncertainties on m; (k) EVPAs (deg); (l) uncertainties on the EVPAs; (m) a and b are the FWHM of the major and minor axes of the modelfit Gaussian components (mas).
ridge-line analysis the HEADACHE 5 python package developed by J. Liu was used.

Images
In Fig. 1

Light curves
In Fig. 2 we show the core region total (upper image) and polarised (lower image) intensity light curves. In each image results at 15, 24, and 43 GHz are reported in single frames from top to the bottom, respectively. We monitored PG 1553+113 between February 2015 to June 2017 with a cadence of about 2 months, and with an approximately six-month gap between the second and third observing runs. Flux densities are reported in Table 2. We find two maxima As shown in the lower image in Fig. 2 we detect polarised emission from the source core region with an average value of 2.6, 2.7, and 2.5 mJy at 15, 24, and 43 GHz, respectively. The polarised emission is also found to be variable, with F var 0.40±0.05 at all three observing frequencies. During MJD 57276 a peak in the polarised emission of 5.6 ± 0.7, 5.5 ± 0.6, and 4.8±0.6 mJy is detected at 15, 24, and 43 GHz, respectively. The fractional polarisation varies in the range 0.7 − 3.1% at 15 GHz, 0.6 − 3.5% at 24 GHz, and 1.1 − 3.7% at 43 GHz. The average fractional polarisation is 1.7%, 2.0%, and 2.2% at 15, 24, and 43 GHz, respectively.
To asses the significance of the variability observed both in the total and polarised intensity light curves, we performed a χ 2 analysis by testing the null hypothesis of non-variability (with the flux density being constant around the average value). In all cases the null hypothesis can be rejected with a confidence level above 3σ.

Jet total intensity ridge line
It is apparent upon visual inspection of the source images ( Fig. 1) that the jet direction varies across epochs. Given that the extended emission of the source on VLBI scales is not prominent and well-defined, it is not straightforward to represent the jet brightness distribution by means of Gaussian model-fit components. In order to constrain the jet direction and to further quantify its variation with time, we calculate the total intensity ridge line in each epoch at 15 GHz after convolving the images with a circular beam with a 0.6 mas radius. We use the method proposed by Pushkarev et al. (2017) by taking azimuthal slices along the jet in a polar coordinate system centred on the core. We look for the half-point that divides the intensity integrated along the arc into two equal sections. The slices along the jet direction are taken in the image plane in stpdf of 0.03 mas and only pixels with a signal-to-noise ratio of S/N > 10 are taken into account. This approach is particularly suitable in the case of a jet showing a limb-brightened structure in total intensity, that is a distinctive feature of HSP blazars (e.g. Giroletti et al. 2008;Piner et al. 2010;Lico et al. 2014). As an example, in the left panel of Fig. 3 we present the 15 GHz total intensity image obtained during MJD 57276, convolved with a 0.6 mas × 0.6 mas circular beam, with the ridge line (yellow line) overlaid. The limb brightened structure is clearly visible. The same image is shown in polar coordinates in the right panel of Fig. 3, where the total intensity limb-brightened structure is also clearly visible (blue colour).
In the upper panel of Fig. 4 we present the ridge lines obtained for all of the 14 observing epochs in a single plot with different colours, where the jet wandering can clearly be seen. Overall, the ridge line extends from 0.9 to ∼ 1 mas from the core region. The jet direction in each epoch is determined by averaging all the angles of each ridge-line point between 0.6 mas from the core region (i.e. beyond the beam radius) and the outermost ridge-line point ( 1 mas). All angles thus obtained are shown in the lower panel of Fig. 4 and reported in Table 3. In principle we could also extract information about the ridge-line curvature. However, given that the jet of PG 1553+113 between 15 and 43 GHz is not particularly extended and is not well defined (core-dominated source), any estimates regarding the curvature would be largely uncertain.
To estimate the angle of the funnel in which the jet wobbles, we first align the 15 GHz images (convolved with a common beam with FWHM 0.6 mas × 1.2 mas) according to the position of the fitted core component, and after subtracting the core component emission we stack all of the residual images. Only those pixels with a S/N > 8 are taken into account. In the residual stacked image (Fig. 5), the jet-emitting region extends up to ∼ 1 mas (∼ 5.2 pc) in a funnel with an angle of φ ∼ 100 deg. A total intensity limb-brightened structure is clearly visible, with the southeast limb being brighter (peak flux density ∼ 3.7 mJy/beam) than the northern one (peak flux density ∼ 2.5 mJy/beam).

Intrinsic polarisation angle and rotation measure
Due to the lack of polarisation calibrators on VLBI scales, in order to obtain the absolute orientation of the EVPAs, a comparison with quasi-simultaneous single-dish or JVLA observations is required. The final EVPA values for the four epochs for which close in time JVLA observations were available (see Sect. 2) are shown in the bottom panel of Fig. 6 for 15 (black circles), 24 (red triangles), and 43 (blue stars) GHz. The EVPAs are found to be variable across the different observing epochs, and change from being roughly aligned (∼ 50 • ) to approximately transverse (∼ 150 • ) to the jet axis.
Because of the effect of Faraday rotation, which occurs when a polarised wave propagates through a magnetised plasma, the observed polarisation angle (χ obs ), at a given observing wavelength λ is rotated with respect to the intrinsic angle (χ int ). Taking this effect into account is essential for obtaining the intrinsic polarisation angle as well as the intrinsic magnetic field orientation. For a Faraday rotating medium external to the emitting region χ obs = χ int + RM × λ 2 , with RM being the rotation measure defined as: where n e is the electron density (cm −3 ), B is the component of the magnetic field parallel to our line of sight (mG), and dl the path length (pc). It is evident, within this theoretical framework, that a linear dependence exists between χ obs and λ 2 . We therefore compute linear fits to our observed EVPAs at 15, 24, and 43 GHz to obtain estimates of χ int and RM. The values of RM and χ int obtained for these four observing epochs are presented in Fig. 6 (top and middle panels, respectively). The average RM value is ∼ −1.0 ± 0.4 × 10 4 rad m −2 , varying between −1.3 and −0.8 × 10 4 rad m −2 . The intrinsic polarisation angle is also variable across the different observing epochs, varying between 135 • ± 7 • (i.e. roughly transverse to the jet axis) and 210 • ± 7 • (i.e. roughly parallel to the jet axis).

Brightness temperature measurements
By fitting the core brightness distribution in the uv-plane with an elliptical Gaussian component, via the modelfit procedure in DIFMAP, we can determine the observed rest-frame core brightness temperature T obs B,vlbi (e.g. Tingay et al. 1998): where S core corresponds to the fitted core flux density (Jy) for a given observing frequency ν (GHz), θ ma j and θ min are the FWHM (mas) of the major and minor axes of the elliptical Gaussian core component, and z is the redshift. The resulting T obs B,vlbi values are reported in Table 5, where T obs B,vlbi (max), T obs B,vlbi (min), and T obs B,vlbi (avg) represent the maximum, the minimum, and average values at each frequency across the different observing epochs.
We also estimate the rest-frame core variability brightness temperature T obs B,var using a variability timescale of τ ∼ 300 days, which is the average time range between the maximum and minimum values observed in the core total intensity light curve: with ∆S max being the difference between the maximum and the minimum core flux density values and d L (m) the luminosity distance (e.g. Liodakis et al. 2017). The T obs B,var values at the three observing frequencies are reported in Table 5 (fourth line).
Since blazar jets are closely aligned to our line of sight, the effects of Doppler boosting must be taken into account when investigating the intrinsic source properties (Urry & Padovani 1995). T obs B,vlbi and T obs B,var are related to the intrinsic brightness temperature T int B via the following equations: where δ = [γ(1 − β cos θ)] −1 is the Doppler factor (with θ being the viewing angle between the jet and our line of sight, β the jet speed in units of the speed of light, and γ = (1 − β 2 ) −1/2 the bulk Lorentz factor).

Discussion
In recent years, several multi-frequency observing campaigns have been devoted to the study of the quasi-periodic variations (on timescales of ∼ 2 years) detected in the γ-ray light curves of the PG 1553+113 (Ackermann et al. 2015). While optical periodic variations on different timescales have been extensively investigated in blazars (e.g. Valtonen et al. 2006;Li et al. 2009;Britzen et al. 2018), at γ-ray energies this has been possible only with the advent of the Fermi-LAT in 2008 thanks to the continuous and systematic sky monitoring in the MeV-GeV energy range. A similar approximately two-year quasi-periodicity (with significance above > 3σ) has also been claimed in a few other blazars such as PKS 0537-441 (   . We emphasise that caution is required when claiming quasi-periodicities on timescales of a few years in blazar γ-ray emission. This is because of the limited γ-ray light-curve duration (based on ∼ 10 years of LAT operational activity), and the red-noise contamination that can easily mimic such short-period variability patterns (e.g. Covino et al. 2019).
One approach to further substantiate the existence of γ-ray periodicity is the use of comparisons to other wavelengths. In the case of PG 1553+113, the optical light curve shows a variability trend in good agreement with the γ rays, confirming the approximately two-year periodic pattern. Moreover, hints of possible periodic behaviour are also found at X-rays, although this is not statistically significant due to the very sparse temporal sampling (Ackermann et al. 2015;Tavani et al. 2018). In this work we extend the investigation to the radio emission from the innermost jet region by means of a dedicated multi-frequency VLBA observing campaign. It is apparent from the 15, 24, and 43 GHz VLBA light curves shown in Fig. 2 that two periods of enhanced radio emission occur around October 2015 (MJD 57312) and June 2017 (MJD 57916), with a minimum appearing around June 2016 (MJD 57558). We note that the spectral index flattens during the periods of enhanced activity, with α 15−43 ∼ 0.2±0.1, while it is steeper during the low-activity state, with α 15−43 ∼ 0.5 ± 0.1 (see Table 3). No clear repeating patterns are discernible in the light curve during the period 2015-2017. Furthermore, the detection of multiple cycles on a longer observing period would be necessary to claim quasi-periodicity in the parsec-scale radio emission.
We note that the 15 GHz VLBA flux density trend over the entire observing period is in good agreement with the light curve obtained from the OVRO monitoring (red triangles overlaid in the top panel of Fig. 2). Indeed, we find F var ∼ 0.13 ± 0.03 for both light curves and the Pearson correlation coefficient between our 15 GHz VLBA and the closest OVRO observations (in time) is r = 0.90, implying a correlation above a 3σ level. These findings allow us to argue that: (1) the overall variability in the radio emission is driven by changes occurring within the VLBI core region, and (2) that the OVRO light curve is representative of the core region flux density modulation. Based on these assumptions, it would appear that despite the presence of rapid variability in the OVRO light curve since 2008 (lower frame in Fig. 7) there is no clear and well-defined periodicity in the radio emission (in contrast to the well-defined periodicity in the γ-ray emission). This is confirmed by applying the weighted wavelet ztransform (WWZ) method (Foster 1996), which is widely used for testing and quantifying quasi-periodic oscillations in blazar light curves (e.g. Hovatta et al. 2008;Ackermann et al. 2015;Bhatta 2017;Tavani et al. 2018;Ait Benkhali et al. 2019). We adopt a WWZ period search window in the frequency range 5.5 − 20.0 × 10 −4 day −1 , with a frequency step of 3 × 10 −5 day −1 , resulting in approximately 50 test frequencies. We choose the minimum frequency in such a way that at least two complete cycles could be detected throuth the full observing time range that covers ∼ 10 yr, while the maximum frequency corresponds to a timescale of 500 days (i.e. about 80 times the mean time separation). The wavelet window width is defined by using a decay constant c = 0.0125. In Fig. 7 we show the OVRO light curve WWZ power (colour scale in the top-right corner) as a function of time (x − axis) and period (y − axis) over the MJD range 54696-58153. From this analysis, no variability pattern is seen in the OVRO light curve over the entire time range, in contrast to that seen in the γ-ray light curve (Ackermann et al. 2015;Tavani et al. 2018). A possible quasi-periodic modulation is found over a limited MJD range (57300-58153) with a period of ∼ 930 days (∼ 2.5 years).
We assess the statistical significance of the WWZ variability timescale by means of Monte Carlo simulations, following the approach proposed by Emmanoulopoulos et al. (2013, and references therein), which takes into account the so-called 'red noise leakage' and aliasing effects due to the sampling properties of the real data sets (i.e. finite length and finite time resolution). Notes. We first generate 2000 light curves with the same power spectrum density (PSD) and power density function (PDF) as the actual observations. We then calculate the WWZ periodograms for the simulated light curves using the same parameters as for the OVRO 15 GHz observations. The averaged WWZ is then compared with the one obtained from the observed light curve. The significance intervals are over-plotted as white contours in Fig. 7, in which the dotted, dashed, and solid lines represent 1σ, 2σ, and 3σ, respectively. These results confirm the lack of a clear periodic variability or pattern that is stable over time at radio frequencies.

A wobbling jet
In general, quasi-periodic MWL emission modulation in blazars has one of two possible origins: (i) The first is a geometrical origin with periodic variations in the Doppler beaming factor produced by periodic changes of the angle between the jet axis and our line of sight. These changes in orientation may be related to the jet precessing or rotational motion, and/or helical structure within relativistic jets (e.g. Camenzind & Krockenberger 1992;Abraham 2000;Stirling et al. 2003;Nakamura & Meier 2004;Rieger 2004;Raiteri et al. 2015). (ii) The second is an intrinsic origin, with quasi-periodic plasma injection into the relativistic jet that is due to quasi-periodic instabilities in the accretion flow, producing the observed variability patterns in the emitted radiation (e.g. Honma et al. 1992;Tchekhovskoy et al. 2011;Pihajoki et al. 2013). In the case of a geometrical effect, the quasi-periodic variability is only observed in the frame of the observer, while in the case of pulsational instabilities in the accretion flow, such quasi-periodic emission variations are also present in the jet comoving frame. In both scenarios the presence of a binary SMBH system is often invoked (e.g. Cavaliere et al. 2017;Sobacchi et al. 2017;Tavani et al. 2018;Caproni et al. 2017;Britzen et al. 2018). As suggested by Begelman et al. (1980), binary SMBHs may produce wiggles, helical motion, and periodic variability in the jets of AGNs as a kinematical consequence of their orbital motion and of jet precession. In order to test the jet precession scenario in PG 1553+113, we constrain the jet position angle across epochs by means of the total intensity ridge line, as described in Sect. 3.3. The jet position angle, as shown in Fig. 4, is found to be variable across the different epochs with a range of values between ∼ 40 • and ∼ 60 • . We assessed the significance of the PA variability by means of a χ 2 analysis, and the null hypothesis of a constant PA equal to the average value of 50 • could be rejected with a confidence level above 3σ. A direct signature of the PA variation across epochs is the large angle (∼ 100 • ) of the funnel in which the jets wander. This was determined via stacking the images of the 15 GHz residuals for all the observing epochs (Fig. 5). The wobbling jet motion indicates that geometrical effects can play a role in the observed emission variability through Doppler boosting modulation. However, we do not find any clear connection between the jet position angle variation and the total intensity emission or with the approximately two-year γ-ray quasi-periodic variability pattern. For this reason, under the assumption that γ rays and radio emission are produced in the same region, we argue that the Doppler boosting modulation alone cannot account for the observed recurrent oscillations in the γ-ray light curve. Additional physical mechanisms are required. Within the framework of a binary SMBH system, one plausible scenario is that the secondary black hole crosses and perturbs the accretion disc of the primary black hole, inducing a temporary enhancement of the accretion rate, which in turn leads to increased jet emission (Sillanpaa et al. 1988;Valtaoja et al. 2000;Caproni et al. 2017;Britzen et al. 2018). Another mechanism that could be responsible for the jet wobbling is the possible shuttle of the core closer to the jet apex, as reported in several studies (e.g. Kovalev et al. 2008;Niinuma et al. 2015). This core-shuttle effect can result from either changes in the physical conditions at the jet base or the ejection of a new component blended within the unresolved core region (e.g. Hodgson et al. 2017;Lisakov et al. 2017).
One possible explanation for the presence of quasiperiodicity in γ rays but not in radio could come from distinct emission regions between the two frequencies. Within this con- text, an additional scenario for producing enhanced γ-ray emission is described by Bosch-Ramon et al. (2012) by means of the possible dynamical interaction of a relativistic jet with matter clumps or the atmosphere of an evolved star in the near vicinity of a SMBH. As the authors of this work point out, in general a considerably amount of dust, stars, and gas tend to cluster in the central region of galaxies and the relativistic jets in AGNs are expected to frequently interact with such ambient material.

Polarisation properties
The polarised emission from the core region (lower image in Fig. 2) is on average more variable (F var ∼ 0.40 ± 0.05) than the total intensity emission (F var ∼ 0.20 ± 0.03), and no variability patterns can be easily recognised during this 2015-2017 VLBA observing campaign. The most notable features in the polarised emission light curves are represented by a ∼ 5.0 ± 0.6 mJy peak reached at all of the observing frequencies during the third observing epoch (MJD 57276), and a second one of ∼ 3.5±0.4 mJy during the fourteenth observing epoch (MJD 57916). Both peaks in the polarised emission are detected during the two periods in which the source is undergoing enhanced activity in the total intensity emission. This behaviour may indicate that the processes producing the increased total intensity emission activity are also responsible for the enhanced polarised emission. The Faraday corrected EVPAs at 15, 24, and 43 GHz, as well as the intrinsic polarisation angle, obtained from the EVPAs versus λ 2 linear fits during MJD 57055, 57078, 57650, and 57755 have been found to vary between being roughly parallel to the jet axis to being roughly transverse, with no clear connection with the total intensity emission modulation. This indicates that the magnetic field configuration is also variable with time. Such behaviour, together with the low degree of polarisation in the core region (on average ∼ 2.5 mJy), could be explained by the presence of multiple polarised subcomponents blended within the beam. The net observed polarisation properties integrated over the VLBI core region result from the sum of the relative contributions of each subcomponent, whose properties vary with time. A&A proofs: manuscript no. Lico_PG1553 The linear dependence found between the observed EVPAs with λ 2 indicate that the Faraday screen is mostly external to the emitting region (e.g. Burn 1966). The RM has been found to vary with time in a range between −1.3 and −0.8 × 10 4 rad m −2 . Given that the jet position angle wobbles across the different observing epochs, we can argue that the RM variations are produced by variations of the path length (d l ) through the external Faraday screen and/or variations in the electron density distribution n e (Stirling et al. 2003). The persistent negative sign in the observed RM across the different epochs, within the scenario in which the Faraday screen is mostly represented by the hot accretion flow or wind outflow, could be explained by a misalignment of the jet axis with respect to the wind axis (Park et al. 2019). If compared to Mrk 421 (e.g. Lico et al. 2014) and Mrk 501 (e.g. Hovatta et al. 2012), the RM value found for PG 1553+113 is the highest measured in a blazar of HSP type to date. However, sample sizes are still not large enough for a statistically significant characterisation of the average RM properties of the HSP blazar population. In this regard, the results of recent and ongoing MWL polarisation VLBA observations will be presented in a further publication.

Brightness temperature
By investigating the brightness temperature properties we can obtain important information regarding the physical conditions within the jet, namely the energy balance within the emitting region. We find that the core brightness temperature decreases as the frequency increases, with T obs B,vlbi (15 GHz) > T obs B,vlbi (24 GHz) > T obs B,vlbi (43 GHz). According to theoretical models (e.g. Marscher 1995), there is a critical frequency ν c below (above) which the brightness temperature increases (decreases) with frequency. Within this scenario, and according to our findings, ν c should be < 15 GHz. Furthermore, by investigating the brightness temperature in a large sample of radio jets in a range of frequencies between 2 and 86 GHz, Lee (2014) found on average ν c ≈ 9 GHz. This decreasing trend of T obs B,vlbi with frequency could indicate that within the two emitting regions probed by the 15 and 43 GHz VLBA observations there is an acceleration of the jet flow with increasing distance from the central engine (see also Melia & Konigl 1989;Marscher 1995;Lee et al. 2016;Nair et al. 2019).
The variability brightness temperature T obs B,var (fourth line in Table 5) for each observing frequency is higher than T obs B,vlbi , as expected because of the different dependence on δ. Given the limitations of the current data sets, we note that for estimating T obs B,var the variability timescale is assumed to be of the order of the average time range between the maximum and minimum values observed in the core total intensity light curve. Even though this approximation is not particularly accurate, it is suitable for the aims of the current analysis. Using Eqs. 4 and 5, with T obs B,vlbi at the peak flux density, we can estimate both T int B and δ: .
We obtain an average value for δ of ∼ 1.4, which is a typical value for HSP objects as obtained from VLBI observations (e.g. Giroletti et al. 2004;Lico et al. 2012;Piner & Edwards 2018). Such a low value for the Doppler factor further supports our conclusion that the Doppler boosting modulation alone cannot account for the observed variability.
The resulting T int B is of the order of 1.5 × 10 10 K, well below the ∼ 5 × 10 11 − 10 12 K inverse Compton limit (Kellermann & Pauliny-Toth 1969), and in agreement with the typical T obs B found in HSP blazars (e.g. Piner et al. 2010;Piner & Edwards 2014;Lico et al. 2016). The physical mechanism generally invoked for explaining the T int B values in HSPs is energy equipartition between particles and the magnetic field. Invoking equipartition yields an upper limit of ∼ 5 × 10 10 K (Readhead 1994). The T int B estimated from our 15, 24, and 43 GHz observations is on average below the equipartition limit, indicating that we are probing an emitting region where the magnetic field energy density (u B ) is higher than the non-thermal particle energy density (u p ). Using equation (5d) from Readhead (1994) we find log(u p /u B ) ∼ −4.6. A similar physical scenario is found in the innermost regions of the radio galaxy M 87, as reported in Kim et al. (2018).

Summary and conclusions
The HSP blazar PG 1553+113 has been observed intensely since an approximately two-year quasi-periodic variability pattern was recognised in its γ-ray light curve (Ackermann et al. 2015). In this work we explored the parsec-scale radio properties of the source by means of a 15, 24, and 43 GHz VLBA observing campaign in total and polarised intensity during the period 2015-2017.
Two periods of enhanced activity around October 2015 and June 2017 emerge from the total intensity light curve with a minimum around June 2016. However, in contrast to the γ-ray emission, no hints of quasi-periodic variability are found in the VLBI emission or in the 15 GHz OVRO light curve over a period covering about 10 years.
We detect polarised emission in the core region (with a polarisation percentage varying in the range ∼ 1 − 4%) that is variable across epochs with no clear recurrent patterns in the light curve. The polarisation angle has been found to be variable across epochs as well, but without any clear connection to the total intensity or polarised emission. We also find a variable RM in the core region, ranging between −1.3 and −0.8 × 10 4 rad m −2 .
The core brightness temperature is found to decrease with increasing frequency, in agreement with Lee et al. (2016), likely suggesting that an acceleration of the jet flow is occurring within the emitting regions probed by the 15 and 43 GHz VLBA observations (Melia & Konigl 1989;Marscher 1995). Using both T obs B,vlbi and T obs B,var we estimate a Doppler factor of ∼ 1.4 and T int B ∼ 1.5 × 10 10 K, indicating that within the emitting region the total energy is dominated by the magnetic field.
By means of a total intensity ridge-line analysis we constrain the jet position angle across the different observing epochs. We find that the jet direction varies in range between ∼ 40 • and ∼ 60 • , indicating that geometric effects could play a role in the observed emission variability through Doppler boosting modulation. However, there is no direct and clear connection between the jet wobbling motion and either the radio flux density or the γ-ray variability pattern, and additional physical mechanisms are invoked. One possible mechanism responsible for the jet wobbling is the core-shuttle effect described by Hodgson et al. (2017) and Lisakov et al. (2017). Another possibility is the presence of a binary SMBH system at subparsec separation, inducing a precessing motion in the jet as well as perturbations in the accretion disc with a consequent modulation in the accretion rate (e.g. Ackermann et al. 2015;Caproni et al. 2017;Tavani et al. 2018). While evidence of binary SMBHs systems at kiloparsec scales have been found in a few tens of objects through direct detection of the two centres, parsec or subparsec systems, which are expected to form quickly according to evolution models, are more elusive. The detection of parsec-or subparsec-scale binary SMBH systems would also be an important assessment of the prediction of coalescence of SMBHs and of consequent gravitational wave production (Begelman et al. 1980;Bhatta 2018;Ait Benkhali et al. 2019). This goal could be achieved in the next few years thanks to the efforts of the Event Horizon Telescope project (Event Horizon Telescope Collaboration et al. 2019), which could potentially allow us to spatially resolve sub-parsec binary SMBH systems.
With the present VLBA monitoring we provide new and valuable information for MWL studies of this peculiar object, furthering our understanding of the physical mechanisms that produce the observed periodicity in the γ-ray emission. This work is intended to be part of a wider and extensive MWL observing program with regular monitoring of the source since 2015 at different frequencies, the results of which will be presented in a dedicated paper (MAGIC Coll. et al., in prep).