Evolution of magnetic field of the Quasar 1604+159 at pc scale

We have analyzed the total intensity, spectral index, linear polarization, and RM distributions at pc scale for the quasar 1604+159. The source was observed in 2002 and 2020 with the VLBA. Combining the MOJAVE results, we studied the evolution of the magnetic field. We detected a core-jet structure. The jet extends to a distance of ~25 mas. The jet shape varies slightly with time. We divided the source structure into the central region and the jet region. In the jet region, we find the polarized emission varies with time. The flatter spectral index values and EVPA direction indicate the possible existence of shocks, contributing to the variation. In the central region, the derived core shift index k_r values indicate that the core in 2002 is close to the equipartition case while deviating from it in 2020. The measured magnetic field strength in 2020 is two orders of magnitude lower than that in 2002. We detected transverse RM gradients, evidence of a helical magnetic field, in the core. At 15 GHz, in the place close to the jet base, the polarization direction changes significantly with time from perpendicular to parallel to the jet direction. The evolution of RM and magnetic field structure are potential reasons for the observed polarization change. The core |RM| in 2020 increases with frequency following a power law with index a = 2.7, suggesting a fast electron density fall-off in the medium with distance from the jet base.


INTRODUCTION
Synchrotron emission is believed to be the main radio emission mechanism of jets of active galactic nuclei (AGN).It results from the acceleration of the ionized electrons in the magnetic (B) field and is highly polarized with theoretical maximum fractional linear polarization of 75 percent (Pacholczyk 1970).
The observed polarized emission strongly connects with the adjacent B field, which is one of the key factors contributing to the propagation of jets.B field wound up by rotating accretion disk or central supermassive black hole (SMBH) extends with jets to distance of parsec (pc) scale even up to kilo-parsec (kpc) scale (Gabuzda et al. 2015a).
Thus, the study of the B field at pc scale is important for the investigation of the physical origin of the propagation of jets.
B field could be detected by linear polarization observation (e.g.Attridge et al. 1999;Gabuzda et al. 2000;Piner et al. 2010;Gabuzda et al. 2014).According to Pacholczyk (1970), the projected components of B field onto the plane of the sky are orthogonal to the polarization angle χ in optically thin regions, while parallel in optically thick regions.This parallel relation happens at an optical depth τ ⋍ 6 − 7 (Wardle 2018).
The information on the B field could also be obtained through Rotation Measure (RM).A set of linearly polarized electromagnetic waves emitted by an AGN jet would experience a rotation of polarization direction when moving through ionized materials, called Faraday Rotation, which could be described with RM, (1) (Jackson 1975), where χ obs is the observed polarization angle, χ o is the intrinsic polarization angle, λ is the wavelength, n e is the density of the electrons in the medium, and B • dl is the projected B filed component on the line of sight (LoS) multiplied by the integral element dl .Thus, the projected components of the B field on the LoS could be detected through RM (Asada et al. 2002;Zavala & Taylor 2005;Hovatta et al. 2012).Further, detection of transverse RM gradients across jets provides strong evidence of helical B field (e.g.Gabuzda et al. 2004;Asada et al. 2008Asada et al. , 2010;;Gabuzda et al. 2008;Kharb et al. 2009;Mahmud et al. 2009;Croke et al. 2010;Hovatta et al. 2012;Gabuzda et al. 2015bGabuzda et al. , 2018)).Clausen-Brown et al. (2011) concluded that a helical B field could also produce significant transverse asymmetries of profiles of fractional linear polarization and spectral index.
Blazars, which include Flat-Spectrum Radio Quasar and BL Lac objects , serve as a strong population of AGN due to the condition that their jet axes orientate close to the LoS of the observer.Emission received would experience an enhancement of power due to the Doppler boosting effect, which makes blazars better candidates to study the B field in jets of AGN.
Source 1604+159 (J1607+1551) has been classified as Low-spectral peaked (LSP) Quasar with redshift z = 0.4965 (Adelman-McCarthy et al. 2008).It has a Luminosity Distance of D L = 2799 Mpc and a linear scale of 6.06 pc/mas.It was detected by both γ − ray detector EGRET (Thompson et al. 1995) and Fermi (Abdo et al. 2010;Nolan et al. 2012).It has flat radio spectral and low radio variability (Mingaliev et al. 2014;Sotnikova et al. 2022).It is a core-jet source, from VLBI up to 8.4 GHz VLA-A scales (JVAS) (Augusto et al. 2006).(Lister et al. 2019) analyzed its jet speed at pc scale and found a maximum 71 ± 14 µas/y (2.09 ± 0.41 c) based on 3 moving features.Pushkarev et al. (2023) and Zobnina et al. (2023) have studied its evolution of linear polarization at 15 GHz spanning from 2009 to 2013.They have found that the jet of 1604+159 has a slight S-shape and the polarization direction or electric vector position angle (EVPA) in general orients towards the jet direction from the VLBI core to a distance of about 4 mas from the core.This means stable toroidal components of the B field exist during the propagation of the jet at pc scale.In this work, we extend the frequency spanning to 4.6 ∼ 43.9 GHz and the time spanning to ∼ 18 years with two multi-frequency observations We further investigate the evolution of the B field and surrounding medium for the quasar 1604+159 through distributions of spectral index, linear polarization, and RM.

Observations
We organized two multi-frequencies dual-polarization observations for the quasar 1604+159 with the Very Long Baseline Array (VLBA), which are the project BH096 (2002-Jul-20) and the project BH228A (2020-Nov-20).
Project BH096 has 3 frequencies centered at 5.0, 8.4, and 15.4 GHz.Each frequency has 4 intermediate frequencies (IFs) and 32 MHz bandwidth.The data rate is 128 Mbits s −1 .The on-source time is 56 mins.
Project BH228A has 7 frequencies spanning 3 bands with 4 in the C band, 2 in the Ku band, and 1 in the Q band, with the DDC personality of the RDBE.Each frequency in the C band has 2 IFs and 128 MHz bandwidth.Frequencies in the Ku band Q band have 4 IFs and 256 MHz bandwidth.The data rate is 2 Gbits s −1 .The on-source time is 92 mins.The information on the observations is listed in Table 1.

Data Reduction
The VLBA raw data was correlated using the correlator at the Array Operations Center of National Radio Astronomy Observatory (NRAO) in Socorro with an average time of 2 s.The initial calibration was conducted using the Astronomical Imaging Processing Software (AIPS) package (Greisen 2003) with the standard procedure.Each frequency was calibrated separately.
For BH228A, the gain values of the C band were not attached correctly.We updated GC tables for the C band data with the appropriate gain values after loading the data into the AIPS with the task 'FITLD' following the suggestion provided by NRAO1 .
After inspection of the data, Fort Davis (FD) and Los Alamos (LA) were chosen as the reference antenna for BH096 and BH228A, respectively.The a-prior amplitude calibrations were done using the information for all antennas in the gain curve (GC), system temperature (TY), and weather (WX) tables to correct the atmospheric opacity.The ionospheric delay was corrected via total electron content measurements from global positioning system monitoring.The phase contributions from the antenna parallactic angles were removed before any other phase corrections were applied.The delay search, bandpass calibration, and RL delay were performed with the source 3C 279 for BH096, and J1751+0939 for BH228A.
The instrumental polarization (D-terms) solutions were solved using the task 'LPCAL'.Source J1407+2827 (OQ208) served as D-terms calibrator for BH096.In BH228A, the compact polarized source J1751+0939 with good parallactic angle coverage during the observation was used.
During the calibration of the EVPA of BH096, the compact polarized source J2253+1608 was used for 5.0 and 8.4 GHz, and 3C 279 for 15.4 GHz.The reference EVPAs were obtained from the Very Large Array (VLA) polarization monitoring program2 and Monitoring Of Jets in Active galactic nuclei with VLBA Experiments (MOJAVE) 2cm Survey (Lister et al. 2018).The nearest epoch for the frequencies 5.0 GHz and 8.5 GHz was 2002 Jul 10.Differences between the reference and our integrated EVPAs were 29.8 • and −41.9 • with errors 2.6 • and 2.0 • for 5.0 GHz and 8.4 GHz, respectively.The nearest epoch for the frequency 15.4 GHz was 2002-July-19.The difference between the reference and our integrated EVPAs was −21.9 • with error 3.0 • (Hovatta et al. 2012).
We arranged two EVPA calibrators for BH228A, source 3C 286 and 2200+420.The source 3C 286 was used for the data of 4.6 and 5.1 GHz (Cotton et al. 1997).The data of 15.2 GHz was calibrated using the MOJAVE 2cm Survey results of the source 2200+420 and the nearest epoch was 2020-Nov-29.After the EVPA calibration of 4.6 and 5.1 GHz, we fitted a linear model to EVPA -λ 2 for 2200+420 (as shown in Fig. 1) and interpolated for other frequencies except for 43.9 GHz data to obtain the correct EVPA.This interpolation method of EVPA calibration was also used in other studies including O'Sullivan & Gabuzda (2009) and Kosogorov et al. (2022).The 43.9 GHz data was calibrated with the results of VLBA Boston U3 .The calibration errors are ∼ 3 • .Detailed information on calibrations is listed in Table 1.
Imaging and self-calibration were carried out in the Difmap package (Shepherd 1997).The I, Q, and U distributions were then obtained using the fully self-calibrated visibilities with natural weighting.For the C band data of BH228A, according to the same suggestion above, as the source 1604+159 is strong enough for self-calibration, we first created a clean model of the source excluding the antennas Pie Town (PT) and North Liberty (NL) in Difmap.Then we computed the self-calibration solutions based on that model with the task 'CALIB' in AIPS, and applied the solutions to all antennas with the task 'CLCAL'.
The distributions of the polarization intensity ( p = Q 2 + U 2 ), polarization fraction (i.e.m = p/I), and EVPAs were calculated pixel by pixel.The corresponding uncertainties were obtained using the approach of Hovatta et al. (2012), which takes into account the effect of the D-terms for each pixel.

Total intensity distribution
Fig. 2 and Fig. 3 show the stokes I distributions for the source 1604+159 at pc scale.There is a core-jet structure.The jet propagates to the direction of ∼ 66 • north by east.In 2002, the jet structure extends to a distance of ∼ 15 mas from the core (5.0 GHz in Fig. 2).With higher sensitivity, in 2020, the jet extending to a distance of ∼ 25 mas from the core is detected (4.6 and 5.1 GHz in Fig. 3).The total intensity map at 43.9 GHz (Fig. 3  and a jet feature (labeled as J) located at ∼ 9 mas from the core.To our knowledge, this is the first time the quasar 1604+159 is detected at 44 GHz.The information on the images is listed in Table 2.
In order to probe whether the shape of the jet varies with time, we collected our 15 GHz data and the Mojave 15 GHz data from 2009 ∼ 2012 (Lister et al. 2018) and used the method described in Cohen et al. (2015); Pushkarev et al. (2017) to derive the jet streamline.This method starts from the core and finds the ridge line points at the same intervals (0.6 mas in this study), which were then connected with a cubic spline.At each ridge line point, the integral of the intensity along a circular arc across the jet, which has the core as its center, is equal on the two sides of the arc.As the method may be affected by the restored elliptical beam (Pushkarev et al. 2017), we restore the median of circular beams derived from the naturally weighted elliptical beams (d = √ major axis × minor axis) for all the 15 GHz images.During the calculation, we blanked pixels with a total intensity < 6σ.
The resultant streamlines are plotted in the Fig. 4. We found the jet propagates along the direction between the core and a feature J.We also found the jet feature J exists at all frequencies in the two epochs, the MOJAVE results from 2009 ∼ 2012 (Lister et al. 2018), and the VLBA Imaging and Polarimetry Survey (VIPS) results (Helmboldt et al. 2007).We believe this is a standing feature detected in most of the observations spanning about 18 years from 2002 to 2020.The information on the standing feature at 15 GHz is listed in Table 3.  , 2, 4, 8, 16, 32, 64, 128).
Interestingly, in the content between the core and the ∼ 6 mas place, there is a hint that the jet rotates clockwise with time from south shown by the jet in 2002 (blue solid line) to north in 2020 (pink solid line).The physical origins behind this phenomenon include the existence of a helical magnetic field as shown by Chen & Zhang (2021), jet procession like M87 (Cui et al. 2023), interaction with the medium surrounding the jet.To note, further long-term observations are needed to confirm whether the rotation does exist and to probe the physical origin.
Although the jet shape slightly changes with time, it goes through the region where standing feature J is located.Besides, maps at the C band in 2020 (Fig. 3) show the eastern jet of 1604+159 continues to propagate through the feature without changing its direction.Thus, we took a straight line through the core and the standing feature J as its jet ridge line (Fig. 4) for the following spectral and polarization analysis.

Core shift and magnetic field strength
As one of the most important properties of a magnetic field, magnetic field strength at 1 pc distance from the central engine could be estimated from core shift measure using our multi-frequency data.According to Blandford & Königl (1979), the core is located at a region where the optical depth, τ , roughly is one.Due to synchrotron selfabsorption, the core at different frequencies is observed in different regions.Assuming a freely expanding, supersonic,  , 2, 4, 8, 16, 32, 64, 128).
narrow, conical jet, the distance of the core from the central engine follows a power law, r core ∝ ν −1/kr obs , where k r = 1 corresponds to the energy equipartition case.
To measure the core shift, we fitted circular or elliptical Gaussian components in Difmap and then obtained the fitted information for components to be used for both epochs, separately.
For data in 2002, feature J and the component (labeled as K as shown in Fig. 2) roughly at 3.7 mas from the core were used.For data in 2020, the feature J is the only component existing at all 7 frequencies and was then used.
The reason for using the combination of circular and elliptical Gaussian components for the model fitting is that it behaves better than using all the circular Gaussian components, and the Root Mean Square (RMS) of the residual maps after fitting is closer to that after clean.Table 4 lists an example of the fitting parameters of components for 4.6 GHz data in 2020.It also provides information (e.g.RMS, χ 2 ) for comparison between using the combination of circular and elliptical Gaussian components for the model fitting and using all the circular Gaussian components., 2, 4, 8, 16, 32, 64, 128).The jet propagates along the direction between the core and the standing feature J.There is a hint for the jet from the core to the ∼ 6 mas place that the jet rotates clockwise with time.2) component used to derive core shift, (3) frequency pair of core shift, (4) magnitude of core shift, (5) direction of core shift, (6) averaged direction of core shift, (7) core shift index, (8) core shift measure, (9) magnetic field strength assuming equipartition, (10) magnetic field strength without equipartition assumption.

Epoch
Component Frequency pair ∆rcore,   The uncertainties of the fitting parameters were obtained following the approach of Lee et al. (2008), which takes into account the effect of signal-to-noise (S/N).For data at the C band, extended diffuse emissions beyond the standing feature J were cleaned to decrease the RMS of the model images.
Core shift vector (|CS| = ∆r core, ν1ν2 ) was measured using the same method of Pushkarev et al. (2012) for frequency pairs including the highest frequency and each lower frequency data for the two observations, separately.This method derives the core shift from CS = IS -OS, where IS is the difference of the image shift (difference of compact component coordinates of images for a given frequency pair), and OS is the difference of the core coordinates.The corresponding uncertainties were derived by using the method described in the study of Chamani et al. (2021Chamani et al. ( , 2023)).In our study, the uncertainty of the core shift is estimated as where σ IS, ν1ν2 , which equals σ 2 , is the uncertainty of aligning images arising from the position uncertainties of the used compact component between frequency pair ν 1 and ν 2 , and σ OS, ν1ν2 is the sum (in quadrature) of the position uncertainties of the core, σ 2 rcore, ν 1 + σ 2 rcore, ν 2 .The resultant core shift information is listed in Table 5.In 2002, the core shift derived from components J and K show similar magnitude but different directions.The directions of core shift derived from K (∼ 15 • of separation from the jet direction) are closer to the jet direction than that from J (∼ 129 • of separation from the jet direction).We noted that the emission of component K is stronger than that of component J, resulting in smaller uncertainties of position.Also, the core shift was derived using vector computing and is small compared with the distances of the two components from the image centers.Thus, the derived core shift (especially its direction) is more easily affected by the larger uncertainty of position and further distance of the component used.With the information on core shift vectors, we obtained absolute projected values and uncertainties (on the mean direction of the core shift vectors) of the core shift and then fitted these values with a power law curve, where ∆r core is the core shift in mas, ν is the observed frequency in GHz, ν ref is the highest observed frequency for a given observation (15.4 GHz for the observation in 2002 and 43.9 GHz for the observation in 2020), a is a fitting coefficient, and k r is the power law index of core shift.The fitted results are plotted in Fig. 5. k r in 2002 has large uncertainties, compared with that in 2020.This may be due to the absence of more samples from other frequencies especially lower frequencies (e.g.L-band).As studied by Sokolovsky et al. (2011), core shift effect is most obvious at lower frequencies and the uncertainty of k r would be significantly decreased by adding L-band data.In 2002, as discussed above, due to the stronger emission and smaller distance to the core of the component K, the k r derived would be more convincing.In 2002, k r is close to 1, the equipartition case.However, k r measured in 2020 shows a larger deviation.The value of k r (lower than 1) and its change with time is very similar to the blazar 3C 454.3 studied by Chamani et al. (2023) with wide multi-frequency and multi-epoch observations, who found the value of k r for 3C 454.3 varies between ∼ 1, and values lower than 1.
Following Lobanov (1998); Pushkarev et al. (2012); Zamaninasab et al. (2014); Zdziarski et al. (2015), we delivered the magnetic field strength at 1 pc using both equipartition and non-equipartition methods for the two epochs.Both assumptions need information on the half-opening angle φ, the Doppler factor δ, and the viewing angle θ.To calculate θ, the Doppler factor δ has to be known.For the source 1604+159, common estimations of δ including from combining X-ray observations with Very Long Baseline Interferometer (VLBI) component fluxes (Ghisellini et al. 1993) and brightness temperature of the source (T b, obs ) (Readhead 1994) need VLBI frequency observed at the turnover frequency and equipartition assumption, which is not the case for the source 1604+159.In 2002, our VLBI observation did not show the turnover, and in 2020, the derived k r suggests deviations from the equipartition.As the source 1604+159 shows low radio variability (Mingaliev et al. 2014;Sotnikova et al. 2022), the variability method (Jorstad et al. 2005;Hovatta et al. 2009) may also not give good estimation of the Doppler factor δ and thus the viewing angle θ.Therefore, we used the same assumption of Pushkarev et al. (2012) to estimate the magnetic field strength.As suggested by Zamaninasab et al. (2014), when the viewing angle could not be estimated accurately, the jet could be assumed by viewed close to their critical angle, θ ≈ Γ −1 (Pushkarev et al. 2012) for the estimation of magnetic field strength.In that case, δ ≈ Γ and Γ = 1 + β 2 app .We also assumed the half-opening angle φ = 0.13Γ −1 (Pushkarev et al. 2012).
If equipartition holds, the magnetic field strength could be estimated through the following equation (4) (Lobanov 1998), where φ is the half jet opening angle, θ is the viewing angle, δ is the Doppler factor, and Ω rν is the core shift measure defined as where D L is the luminosity distance in parsec.
We found in 2002, with core shift derived from component K, both the values of the B 1 pc measured by using the equipartition and non-equipartition methods are close to each other in magnitude.However, in 2020, the nonequipartition case shows the strength of B 1 pc two orders of magnitude lower than the equipartition case, consistent with the magnetic field study of part of 59 blazars for the two cases by Zdziarski et al. (2015).The B 1 pc derived for 1604+159 also shows change with time.A more comprehensive discussion will be presented in Section 4.

Spectral index
To avoid spurious gradients, images need to be aligned before the construction of spectral index maps due to loss of absolute position after self-calibration for each frequency.Two main methods can be used to align images.One is using bright and compact components in the jet assuming their optically thin nature (e.g.Lobanov 1998;Marr et al. 2001;Kovalev et al. 2008;Sokolovsky et al. 2011;Hovatta et al. 2012).The other is two-dimensional (2D) cross-correlation to obtain its best core-shift for each frequency through calculation of the largest correlation coefficient of optically thin jet for each frequency (Walker et al. 2000;Croke & Gabuzda 2008).As the jet of 1604+159 at pc scale has several bright and compact components detected (Fig. 2 and 3), the former aligning method was adopted.
Source 1604+159 is a LSP Quasar with a spectral peak located at a frequency < 10 14 Hz.Compared with the optically thin nature of the jet region, the spectra in the central region may behave more complexly.To determine the frequency intervals for spectral index maps, we re-produced I distributions for each frequency with the same pixel size, shifted image centers to the centroid of the standing feature J, and converted the images to the beam size of the lowest frequency image for both epochs.Then 5 points (roughly the same interval of ∼ 0.6 mas) were sampled along the jet direction from the central region, shown in Fig. 6.
As seen in Fig. 6, in 2020, the total intensity of 4.6-5.1-6.0GHz and 12.2-15.2-43.9GHz follow a power law for the 5 sampled points in the central region.However, the total intensity at 7.8 GHz shows more complex behavior.It drops earlier than that at 12.2 GHz and higher frequencies (as shown from (1)-( 5) in Fig. 6) when being away from the jet base along the jet direction.Based on that, spectral index maps were constructed for the following frequency intervals: 4.6-5.1-6.0,6.0-7.8,7.8-12.2,and 12.2-15.2-43.9GHz in 2020; 5.0-8.4 and 8.4-15.4GHz in 2002.with 3σ × (−1, 1, 2, 4, 8, 16, 32, 64, 128).Take 4.6-5.1-6.0GHz data in 2020 as an example: The (u,v)-coverage was first flagged to the same range in which the minimum UV distance is the smallest value of the highest frequency and the maximum is the largest one of the lowest frequency.Then, maps were restored to the lowest frequency beam size.Last, the map center of each frequency was shifted to the position of the standing feature J obtained from model fitting to align the images, which was also used by (Kosogorov et al. 2022).Other frequency intervals were treated with the same process.
To derive the spectral index and its uncertainty pixel by pixel, we used the same method of Hovatta et al. ( 2014) .More specifically, the spectral index distribution (S ∝ ν α ) of each frequency interval was obtained by fitting a power-law to the total intensity distribution of each frequency in the interval.When fitting the spectral index, we added in quadrature a calibration uncertainty of 0.1.Pixels were clipped when the total intensity level was less than 3σ ν at the corresponding frequency, where where σ rms is the thermal noise of the image taken at a corner of each map.We also blanked pixels with uncertainty larger than 1.The resultant map is shown in Fig. 7.
To analyze the spectral distribution along the source, we extracted the spectral index values along the ridge line, as shown in Fig. 8.We divided the emission region into two regions based on the polarization distribution in 2002 (Fig. 2): the central region including the core and upstream of the jet, and the jet region.
In the central region, the core in general has flat spectra.For frequency intervals 8.4-15.4GHz in 2002 and 7.8-12.2GHz in 2020, the spectral index experiences a smooth change from a slight inverted value (> 0) to a flat value (∼ 0) to a steep value (< 0).This smooth transition is due to the convolution with the beam, but the change is intrinsic to the source, as shown in the simulations by Hovatta et al. (2014), indicating that the central region of the source 1604+159 experiences a spectral change from synchrotron self-absorbed to flat to optically thin.This change with distance also reflects the opacity change with distance in the core.
For the jet region, Hovatta et al. (2014) measured spectral index values from the jet ridge lines.They found quasars have a mean spectral index of −1.09 ± 0.04.Compared with this mean value, we noted that, in the two epochs, regions with bright, compact jet features have larger spectral values including the component K and the standing feature J. Assuming a homogeneous synchrotron for these regions, the spectral index α = 1−p 2 , where p is the index of the energy distribution of electrons N (E) = N 0 E −p , where E is the electron energy.The increased spectral indices indicate the increase of electron energy in the compact features, which suggests that electrons in the regions are accelerated.This supports the existence of shocks in these flatter regions.The standing feature J detected in most of the observations from 2002 to 2020 provides another evidence of the existence of shock in this region.

Linear polarization distribution
Fig. 2 and Fig. 3 also show linear polarization and fractional linear polarization (m) distributions for 1604+159.The polarization pixel was clipped when its intensity was below 3σ p for each frequency, where σ p is derived according to the polarization uncertainty analysis of Hovatta et al. (2012).
The central region contains polarized emissions from the core and the upstream of the jet, the averaged m is 5.3% at 5.0 GHz, 6.9% at 8.4 GHz, and 10.2% at 15.4 GHz in 2002, indicating the depolarization effect due to reasons including Faraday Rotation and opacity effects.In 2020, the averaged m shows a similar depolarization effect, as shown in Table 2.
In the jet, polarized emissions are centered in areas with bright, compact features.We found the fractional polarization is higher than that in the central region.More and more evidences of increase of fractional polarization with distance from the jet base indicate that magnetic field becomes more ordered along the jet (Cawthorne et al. 1993;Lister & Smith 2000;Lister 2001;Lister & Homan 2005;Pushkarev et al. 2023), which may due to reasons including small Faraday rotation (Hovatta et al. 2012) and spectral ageing effect (Kardashev 1962).
In the core, observations at different frequencies reveal emissions from different regions, with higher frequencies corresponding to regions closer to the jet base.The surrounding medium of the core at different frequencies may also have diverse properties including spectral index and Rotation Measure, which may lead to complex distribution of the observed polarization.To analyze the observed EVPA along the jet, we obtained the distributions of EVPA minus jet direction along the jet ridge line for the two observations, as shown in Fig. 9.
For the central region, Pushkarev et al. (2023) has analyzed the EVPA orientation to the jet direction along the ridge line for blazars.They found that for the core of quasars the |EVPA − jet PA| has a clear predominance of a peak near 90 • and a weaker secondary peak at values close to 0 • , in contrast to BL Lacs which show more predominant EVPAs aligned with the jet direction.As seen in Fig. 9, the core of 1604+159 has a simple consistent polarization direction basically along the jet except for the 15 GHz in 2020.In 2020, a complex variation of the EVPA with distance appears in the 15 GHz core, which experiences a rapid change from perpendicular (∼ 90 • ) to parallel (∼ −25 • ) to the jet direction starting at an area close to the jet base.In contrast, in 2002, the EVPA of the 15 GHz core in general orients toward the jet direction.This suggests that the polarization of the 15 GHz core evolves with time, confirming the results of (Lister et al. 2018) which have shown a similar large change of EVPA with time in the area close to the jet base of 1604+159.
In the jet region, Pushkarev et al. (2023) also found that the EVPAs in quasars tend to be preferentially transverse to the jet in the outer ridge line regions.For the jet of 1604+159, we did not detect the EVPA transverse to the jet but had a large separation of ∼ ±40 • from the jet direction, where positive and negative signs are counterclockwise and clockwise, respectively.Compared with that of (Lister et al. 2018) and Pushkarev et al. (2023), we also found the observed EVPA changes with time in the jet.The origin of the change of EVPA in the central region and the jet region will be discussed in Section 4.

Rotation Measure
Our broad-band frequency (4.6 ∼ 43.9 GHz) data in 2020 allows us to probe RM at different scales from different frequency intervals.To determine the frequency intervals, we imaged polarization distribution for each frequency observed in 2020 with the same pixel size, shifted image centers to the centroid of the standing feature J, and converted the images to the beam size of the 4.6 GHz image.Then a point was sampled from the central region in the images of the seven frequencies, and the EVPA distribution against squared wavelength was obtained for the point as shown in Fig. 10.Observation in 2002 was done in the same process.
As can be seen from Fig. 10, in 2020, the EVPA stays stable at the C band and decreased at 12.2 GHz following a non-linear behavior, with larger EVPA rotation at larger frequencies, leading to a possible relation between RM and frequency, which will be analyzed and discussed in Section 4. In 2002, the EVPA seems to show a more simple rotation.This may be due to the absence of observation at other frequencies or the existence of a simple external medium.Given its frequency range and sample number (5.0, 8.4, and 15.4 GHz), this simple rotation of EVPA could be fitted by a straight line.The foreground RM arising by our Milky Way for the source 1604+159 is negligible because it is located in one of the 'holes' where the amplitude of the median RM is consistently below 1 rad m −2 (Taylor et al. 2009).Thus, the rotation of EVPA shown in Fig. 10 probably reflects the intrinsic physical condition in the central region of 1604+159 and its environments.
We then selected four frequency intervals for the observation in 2020 to construct RM distributions which include 4.  GHz.RM distribution in 2002 was also constructed.
The multi-frequency data allows us to access the Rotation Measure (RM) according to the equation 1. Pixel was blanked when the chi-squared (χ 2 ) exceeded 3.84 (5.99), corresponding to a 95% confidence limit for 1 (2) degree of freedom.The χ 2 is from the standard equation: where O i are the observed values, E i are the expected values of the fitting model, and σ i are the EVPAs variances (e.g.Press et al. 1992;Hogg et al. 2010).
When constructing RM distributions for the frequency intervals 12.2-15.2GHz and 15.2-43.9GHz, we minimized the modulus of RM by shifting the EVPA at one frequency by nπ to solve the EVPA ambiguity problem (Kosogorov et al. 2022).The resulting RM distributions are plotted in Fig. 11, Fig. 12, and Fig. 13.
RM serves as an important tool to probe the B field since it carries information on the projected B field on the LoS, as shown by equation 1.Detection of transverse RM gradients with significance of 3σ and above provides strong evidence for the existence of helical B field (Hovatta et al. 2012;Gabuzda et al. 2018).The majority of transverse RM gradients were found at central regions including cores and upstream of the jets of blazers (e.g.0059+581, 11240059+581, -186, 19080059+581, -201 (Gabuzda et al. 2018))).For the observation in 2002, we detected a transverse RM gradient with 2σ significance at the core region of 1604+159, as shown in Fig. 11.This corresponds to a confidence level > 90% (Hovatta et al. 2012).In 2020, at almost the same place, we detected a transverse RM gradient with 3σ significance within a similar frequency interval, as shown in Fig. 12.The detection suggests a high possibility of the existence of a helical B field in the emission region of the core reflected by the frequency intervals (i.e.5.0 ∼ 15.4 GHz in 2002 and 6.0 ∼ 12.2 GHz in 2020).Although the helical field is believed common in the jet of AGN, only ∼ 50 jets with firm transverse RM gradients are detected.Most jets in AGN have their transverse structure unresolved and not enough polarized emission detected, which makes it difficult to obtain firm transverse RM gradients.This needs a higher resolution and sensitivity array like space VLBI.
When fitting RM, the RM corrected EVPA, which represents the intrinsic polarized direction and thus the projected structure of the B field on the plane of the sky, could also be obtained.Fig. 11,Fig. 12 and Fig. 13 also provide the RM corrected EVPA for the two observations.As can be seen, almost all the derived RM corrected EVPA distributions in the central region are along the jet direction except for that in 12.2-15.2GHz in 2020, which shows a very large rotation of EVPA along the jet direction.The corresponding RM distribution also shows a large change of magnitude with the further distance the smaller the values.This RM and intrinsic EVPA result in the observed polarization in the central region of 15.2 GHz in 2020 (as described in section 3.4.1 and shown in Fig. 9).The origin behind this change may lead to the evolution of the observed polarization at 15 GHz detected by our observations and MOJAVE.

DISCUSSION
In this section, we discuss the evolution of the magnetic field for 1604+159 at pc scale.
4.1.The magnetic field in the Jet MOJAVE monitored the source 1604+159 for 6 epochs from 2009 to 2013 at 15 GHz (Lister et al. 2018) and studied the source with stacking (Pushkarev et al. 2023).Stacking has the advantage of improving the image sensitivity and providing a complete jet cross-section structure as much as possible by utilizing existing data.Through distributions of EVPA and fractional polarization, stacking shows long-term persistent magnetic field configuration.We revisited the stacking fits files derived by Pushkarev et al. (2023) for 1604+159, clipped pixels which have polarization magnitude below 5σ (σ = 5 × 10 −5 Jy beam −1 for this case), and obtained 4 transverse slices of fractional polarization roughly normal to the jet direction in the jet.Fig. 14 plots the resultant slices.Stacking results for this source have shown a stable linear polarization structure extending to ∼ 4 mas from the core for 1604+159 .This polarization gradient basically orientates along the jet direction from the core to the jet, indicating a global toroidal B field component existing in the jet.In addition, the northern side of the jet has higher fractional polarization than that of the southern side, indicating a helical B field rather than a toroidal B field which could produce a symmetrical transverse fractional polarization profile.If a helical B field is in the jet, the side with the direction of B field close to the LoS would have lower fractional polarization than the other side which has the direction of B field close to the plane of the sky (Clausen-Brown et al. 2011).
Higher fractional polarization at one edge of the jet could also be explained by the interaction with ambient media (Gabuzda 2021).However, if a jet propagates through its ambient media, interaction with ambient media would produce a longitudinal B field (i.e.orthogonal linear polarization) (Laing & Bridle 2014), which is not detected in the stable polarized gradient.Thus, the helical B field would be a more natural physical origin of the stable polarization structure.But we could not rule out the interaction possibility.
According to the study of Murphy et al. (2013), given a helical B field, the distribution of intrinsic EVPA depends on the pitch angle of the B field and viewing angle.As the jet of 1604+159 has low RM mediums (Fig. 11, Fig. 12 and Fig. 13) surrounding the jet, the observed EVPA at frequencies above 4 GHz could be used to represent the intrinsic polarized direction.The fact that the jet does not change its shape results in the expectation that the viewing angle would be unchanged between 2002 and 2020.The pitch angle also appears unchanged.Assuming a global helical magnetic field in the jet region, it is expected that a stable longitudinal EVPA distribution would be observed as the jet propagates basically along a line.Comparing with the linear polarization results from 2009 to 2013 (Pushkarev et al. 2023) , we found that the polarization in the jet region changes with time, as can be seen in Fig. 9. Our results show that the polarization in the jet observed in 2002 and 2020 differs from the jet direction with a relatively large separation of ∼ ±40 • or even larger in some places.We noted that the polarization of the jet observed in the two epochs is centered in bright, compact features, which could not be naturally explained as being associated with helical B field (Gabuzda 2021).The large separation of EVPA suggests the helical B field is severely distorted in the regions with compact features.As discussed in section 3.3, shocks may exist in these regions with bright, compact features detected on the two epochs.Shock compressing the B field leads to the B field along the shock front (Laing 1980) .With the relation of 90 • between the EVPA and the B field in an optically thin area (Pacholczyk 1970), the polarization in the jet region in 2002 and 2020 indicates shocks are oblique, consistent with the study of Lister & Homan (2005) who found the shocks could be oblique in quasars.Our results support the point that shocks play an important role in the variability of jet emission, especially polarization variability (Gabuzda 2021).Before exploring the origin of the significant evolution in EVPA, we investigate the potential magnetic field structure in the central region.In section 4.1, we find evidence for the presence of a helical field in the jet region.Here, we provide more observed evidence for the existence of a helical field in the core region.First, the sign changed transverse RM gradients are detected in 2002 and 2020 (Fig. 11 and Fig. 12).Second, the asymmetry of fractional polarization across the jet is reflected in Fig. 14.As pointed out by Murphy et al. (2013), the side when the fractional polarization (Fig. 14) is lower has a larger magnitude of RM (Fig. 11 and Fig. 12), which we do observe in the southern side of the core.Third, the RM corrected EVPA in the core has longitudinal distribution (Fig. 11, Fig. 12, and Fig. 13) except for the 12.2-15.2GHz which will be discussed in the following content.Compared with the results by Murphy et al. (2013), the longitudinal EVPA indicates a large intrinsic pitch angle (∼ 80 • ) in the jet rest frame assuming the polarized emission is dominated by optically thin regions.Therefore, the potential magnetic field structure in the central region of 1604+159 could be a helical field with a large intrinsic pitch angle.
Under external screen assumption, the factors contributing to the evolution of observed EVPA at 15 GHz are RM and intrinsic EVPA (as indicated by equation 1), both of which have a close relation with the B field.
In 2002, the 15 GHz related intrinsic EVPA shows distribution parallel to the jet direction (Fig. 11).However, in 2020, the intrinsic EVPA reflected by 12.2-15.2GHz behaves more complex (Fig. 13).Here, we provide two possible explanations.One is that, as discussed by Gabuzda (2021), the parallel relation between the EVPAs and projected B on the plane of the sky happens at an optical depth τ ⋍ 6 − 7 (Wardle 2018), and the core observed in jet often locates at place having the optical depth ∼ 1 (Blandford & Königl 1979).The observed polarized emission thus is from optically thin regions.Therefore, if this is the case, the change in EVPA reflects a change of the inherent B field in the region of emission at 15 GHz.This complex distribution indicates that the configuration of the helical magnetic field breaks in the region of emission at 15 GHz.
Another is that polarized emission from optically thick regions could not be negligible.For the environments surrounding the jet, 15 GHz related RM in the core in 2020 is much larger than that in 2002 (Fig. 11 and Fig. 13).As the magnetic field strength measured in 2020 is ∼ 2 orders of magnitude smaller than that in 2002, the much larger RM in 2020 results from much higher electron density assuming the same integral region, which leads us to consider the contribution of polarized emission from optically thick regions.A blending of emissions from both optically thin and thick regions results in the structure of the magnetic field being undetermined.However, no matter which scenario is more appropriate, the change of intrinsic EVPA with time at the place close to the jet base reflects the change in the physical environment with time.As the inner jet does not change its direction (Fig. 4 ), this large change (with time) of the B field structure may result from the rapid change of the environments surrounding the jet and/or central engine (e.g. the SMBH or/and the accretion disk) producing the B field.
Under the equipartition assumption, the core RM follows the relation |RM| ∝ ν a , where a is the index in the electron density n e power-law change as a function of distance r from the black hole, n e ∝ r −a (Jorstad et al. 2007).The value of a tells us the potential physics condition in the core.For example, a ∼ 2 is expected for a conically expanding jet having a boundary layer, which is the case of 1156+295, 2007+077, and 3C 273 (O'Sullivan & Gabuzda 2009;Hovatta et al. 2019).
We took the averaged RM and uncertainty in the central region as the core RM and its uncertainty for each frequency interval, listed in Table 6.Then we obtained the index a = 2.7 ± 0.5 from fitting, larger than 2, suggesting the core in 2020 has faster electron density fall-offs in the medium with distance from the jet base.High index values are also observed in the core of 0954+658, 1418+546 (O'Sullivan & Gabuzda 2009), 3C 279, OJ 287, CTA 102, 1749+096, 0235+164, and BL Lac (Park et al. 2018).
The Large value of a indicates some assumptions for a ∼ 2 break, including magnetic field structure and equipartition.Although the value of 2.7 in 2020 does not significantly deviate from 2 if considering 2σ uncertainty, we do find other evidence that some assumptions break.The magnetic field of the core reflected by the RM corrected EVPA distribution for 12.2-15.2GHz in 2020 (Fig. 13) has a more complex structure than that in 2002, a simple helical field.The core shift index k r value of 0.32 suggests the extent of deviation from equipartition is more significant than in 2002.
We have observed significant behavior differences in linear polarization at 15 GHz and in the Faraday Rotation medium in the core in the two epochs.These variations indicate that the B field wound up by the central engine in the core correlates with the medium wrapping up it.Linear polarization at 15 GHz collected spanning ∼ 18 years from 2002 to 2020 suggests the core of 1604+159 experiences a complex evolution of physical condition including the B field and the medium near the jet base.
Other blazars are showing significant evolution of B field and RM.O'Sullivan & Gabuzda (2009) found no good linear λ 2 fits for the source 1749+096 observed in 2006 July, while observation two weeks earlier detected core RM (Hovatta et al. 2012).One possible explanation for this variation is a possible flare in the source affecting the polarization results.Mahmud et al. (2009) found a surprising reversal of Faraday Rotation gradients in the jet of the BL Lac object B1803+784.One of the most likely origins they explained is a 'nested-helix' magnetic field structure, resulting from magnetic-power-type scenarios (Lynden-Bell 1996).The origins of RM change with time are different for these sources.Whether the differences are common in AGN needs studies of the evolution of RM with a large number of sources.

CONCLUSIONS
We have analyzed the total intensity, spectral index, linear polarization, and RM distributions at pc scale for the quasar 1604+159.The source was observed at 5. 0, 8.4, and 15.4 GHz in 2002, and 4.6, 5.1, 6.0, 7.8, 12.2, 15.2, and 43.9 GHz in 2020 with the VLBA.Combining the MOJAVE polarization results at 15 GHz from 2009 to 2013, we studied the evolution of the magnetic field spanning ∼ 18 years.Based on the linear polarization distribution in 2002, we divided the source structure into the central region and the jet region.The conclusions are summarized as follows: 1. We detected a core-jet structure.The jet extends to the distance of ∼ 25 mas from the core at a direction of ∼ 66 • north by east.The shape of the jet derived from 15 GHz data varies slightly with time and could be described by a straight line.
2. The core shift index k r in 2002 indicates the core is close to equipartition, while in 2020 the core shows a large deviation.The measured magnetic field strength in 2020 is two orders of magnitude lower than that in 2002.
3. We find the polarized emission in the jet region varies obviously with time.The polarization has a large direction separation of ∼ ±40 • from the jet direction in the two epochs.These polarized emissions are centered in regions with bright, compact jet features, which have spectral index values along the jet ridge larger than the mean spectral index of quasars in the jet.The flatter spectral index values indicate the possible existence of shocks in these bright, compact jet features, contributing to the variation of polarization in the jet with time.The oblique EVPAs indicate the shocks are oblique.
4. In the central region, the spectral index shows flat spectra and the fractional polarization is low compared with that in the jet region.The polarized emission orientates in general towards the jet direction.At 15 GHz, in the place close to the jet base, the EVPA changes significantly with time from perpendicular to parallel to the jet direction.We find in 2002 the EVPA at 15 GHz in this place is basically along the jet and the rotation of that at the three frequencies could be well fitted to a linear model (external medium screen).In 2020, the EVPA at this place becomes perpendicular, and the rotation of that behaves more complex.
5. We detected transverse RM gradients for the two observations in the core, suggesting the existence of a helical field.The 15 GHz related RM in the core in 2020 shows a large change with distance.The RM corrected EVPA distribution in 12.2-15.2GHz in 2020 reflects a complex magnetic field structure in the core of 15 GHz.The core |RM| increases with frequency following a power law with index a = 2.7 ± 0.5, showing faster electron density fall-offs in the medium with distance from the jet base than a sheath surrounding a conically expanding jet in which the power law index is ∼ 2.

Figure 1 .
Fig.2and Fig.3show the stokes I distributions for the source 1604+159 at pc scale.There is a core-jet structure.The jet propagates to the direction of ∼ 66 • north by east.In 2002, the jet structure extends to a distance of ∼ 15 mas from the core (5.0 GHz in Fig.2).With higher sensitivity, in 2020, the jet extending to a distance of ∼ 25 mas from the core is detected (4.6 and 5.1 GHz in Fig.3).The total intensity map at 43.9 GHz (Fig.3) shows the core

Figure 4 .
Figure 4. Streamline derived to represent the jet shape for 15 GHz data from 2002 to 2020 with 2002 in blue and 2020 in pink.Contours are the 15 GHz total intensity (2002 in blue, 2020 in grey) with 3σ × (1, 2, 4, 8, 16, 32, 64, 128).The jet propagates along the direction between the core and the standing feature J.There is a hint for the jet from the core to the ∼ 6 mas place that the jet rotates clockwise with time.

Figure 5 .
Figure 5. Core shift as a function of frequency with 15.4 GHz and 43.9 GHz as a reference frequency for observations on 2002-Jul-20 and 2020-Nov-20, respectively.In 2002, compact components K and J were used.

Figure 6 .
Figure6.A plot of total intensity against frequency for 5 points (roughly the same interval ∼ 0.6 mas) sampled along the jet direction from the central region (shown in the last map) for the multi-frequency observations in 2002 (blue) and 2020 (orange).The titles (1)-(5) indicate the positions with larger numbers corresponding to the distance further from the jet base.The spectra of the sample show complex change, indicating a potentially complex physical environment in the central region for the source.

Figure 9 .
Figure 9. Distributions of EVPA − jet direction along the jet for observations on 2002-Jul-20 and 2020-Nov-20.The black scatter points are derived from the MOJAVE stacked results at 15 GHz between 2009-2013.The black dashed line divides the source into the central region and the jet region.

Figure 10 .
Figure10.A plot of EVPA vs λ 2 for a point sampled from the central region for the multi-frequency observations in 2002 (blue) and 2020 (orange).In 2002, EVPA changes linearly with λ 2 .In 2020, the behavior of EVPA is more complex.

Figure 11 .
Figure 11.Distributions of RM and RM corrected EVPA for observation on 2002-Jul-20 with a 2σ transverse RM gradient extracted along the transverse slice (blue solid line) in the RM map.The label 'S' is the start point of the transverse RM gradient.Color is the RM and RM error.Contours are the 15.4 GHz total intensity with 3σ × (−1, 1, 2, 4, 8, 16, 32, 64, 128).

Figure 12 .
Figure12.Distributions of RM and RM corrected EVPA at frequency interval 12.2-15.2GHz for observation on 2020-Nov-20 with a 3σ transverse RM gradient extracted along the transverse slice (blue solid line) in the RM map.The label 'S' is the start point of the transverse RM gradient.Color is the RM and RM error.Contours are the highest frequency total intensity at given frequency intervals with 3σ × (−1,1, 2, 4, 8, 16, 32, 64, 128).

Figure 14 .
Figure 14.MOJAVE stacked image for 1604+159 and transverse slices of fractional polarization.Color is fractional polarization.

4. 2 .
The magnetic field in the central region While the central region at 15.4 GHz in 2002 remains EVPAs along the jet direction, it experiences a smooth change from perpendicular to parallel to the jet along the jet direction beginning at a place close to the jet base(Fig.9).This change has also been detected by MOJAVE at epoch 2009-Aug-19 and 2012-Jul-12 4 .The Study of linear polarization variability (Zobnina et al. 2023) has shown that in the place close to the jet base, 1604+159 in 2009 ∼ 2013 has a large change of EVPA with time and the σ EVPA ∼ 50 • .This σ is large in the population of quasars even in the blazars, suggesting the special physical condition of this source.The significant difference in EVPA in the central region between 2002 and 2020 confirms the large σ.

Table 1 .
Observation information

Table 2 .
Image statistics summary.

Table 3 .
Information of the standing feature J at 15 GHz