The NANOGrav 12.5-Year Data Set: Polarimetry and Faraday Rotation Measures from Observations of Millisecond Pulsars with the Green Bank Telescope

In this work, we present polarization profiles for 23 millisecond pulsars observed at 820 MHz and 1500 MHz with the Green Bank Telescope as part of the NANOGrav pulsar timing array. We calibrate the data using Mueller matrix solutions calculated from observations of PSRs B1929+10 and J1022+1001. We discuss the polarization profiles, which can be used to constrain pulsar emission geometry, and present both the first published radio polarization profiles for nine pulsars and the discovery of very low intensity average profile components ("microcomponents") in four pulsars. Using the Faraday rotation measures, we measure for each pulsar and use it to calculate the Galactic magnetic field parallel to the line of sight for different lines of sight through the interstellar medium. We fit for linear and sinusoidal trends in time in the dispersion measure and Galactic magnetic field and detect magnetic field variations with a period of one year in some pulsars, but overall find that the variations in these parameters are more consistent with a stochastic origin.

netic field parallel to the line of sight for different lines of sight through the interstellar medium. We fit for linear and sinusoidal trends in time in the dispersion measure and Galactic magnetic field and detect magnetic field variations with a period of one year in some pulsars, but overall find that the variations in these parameters are more consistent with a stochastic origin.
Keywords: pulsars: general, ISM: magnetic fields, techniques: polarimetric Pulsars are highly-magnetized, rapidly-rotating neutron stars that emit electromagnetic radiation that sweeps across our line of sight as they rotate. In addition to being laboratories for study themselves, pulsars are useful in probing the properties of the interstellar medium (ISM). As the radio waves from a pulsar traverse the Galaxy, they experience Faraday rotation, which is a frequency-dependent rotation of the polarization position angle by the Galactic magnetic field.
where e is the charge of the electron, λ is the wavelength of the radio waves, m e is the mass of the electron, c is the speed of light, n e is the free electron density along a line of sight l, d is the pulsar distance, and B is the estimate of the electron-density-weighted average (Galactic) magnetic field (in cgs units). The degree to which the pulsar's radio waves are rotated is called the rotation measure (RM), where We can also measure directly from radio observations the dispersion measure (DM), which is the integrated free electron density along the line of sight and varies with the observational frequency as ν −2 . We can then calculate the parallel component of the magnetic field along the line of sight using both the DM and RM as where RM is in rad m −2 and DM is in pc cm −3 . When the radio waves reach the receiver on the telescope, the telescope's response alters the components of * NANOGrav Physics Frontiers Center Postdoctoral Fellow the waves. These components can be described by the Stokes vector where Stokes I is the total intensity, Stokes Q and Stokes U form the linear polarization L = Q 2 + U 2 , and Stokes V is the circular polarization intensity. Using the International Astronomical Union (IAU)'s circular polarization sign convention, right-handed circular polarization is positive (corresponding to a clockwise rotation of the position angle) and left-handed circular polarization is negative (corresponding to a counterclockwise rotation of the position angle) (Stokes 1851). The total amount of polarized emission can be described by the latter three Stokes parameters, P = √ L 2 + V 2 . The orientation of the linearly polarized radio waves emanating from the pulsar can be described by the position angle (PA) of the linearly polarized emission: The polarization angle is quoted using the IAU convention with the polarization angle increasing in the counterclockwise direction. We can solve for the telescope's response to the incoming radio waves where S meas is the Stokes vector measured at the telescope, S src is the Stokes vector of the incoming radio waves, and M is the Mueller matrix, which depends on the ellipticity of the receiver arms, non-orthoganality of the receivers, the differential gain, and the differential phase of the receiver (see Heiles et al. 2001 for more details). By observing a strongly-polarized source through a series of telescope azimuth angles, Mueller matrix elements for a given telescope and observing system can be determined. The Mueller matrix can then be used to correct other observational data and recover the intrinsic Stokes parameters of the source under observation. We can determine the Mueller matrix for a certain receiver by taking a long observation of one pulsar and tracking it across the sky. By doing this for multiple epochs, we can judge the stability of the receiver by observing how the solutions change over a long period of time.
Pulsars are highly-polarized sources, and the position angle can vary across the pulse phase. For many pulsars, this follows an S-shaped curve, interpreted through the rotating vector model (RVM; Radhakrishnan & Cooke 1969) as the observer's line of sight traversing a conal emission beam, with radio emission originating from the open magnetic field lines. The position angle is measured with respect to the magnetic axis such that it will rotate through the pulse by at most 180 • .
Numerous polarization studies on millisecond pulsars (MSPs) Xilouris et al. 1998;Stairs et al. 1999) have demonstrated that most MSPs have more complex position angle curves which are notoriously difficult to fit to this model (Craig 2014) (Stairs et al. 1999), and ). This is due to the recycled nature of MSPs, creating a complicated field configuration and a reduction in the magnetic field strength, resulting in much smaller period derivatives than canonical pulsars.
Millisecond pulsars also feature emission over a large portion of the profile, with more complex profiles and less profile evolution with frequency than canonical pulsars . Geometric arguments imply that pulse widths should vary as the inverse square root of the period (Rankin 1993). In addition, the beams of millisecond pulsars are wider than those of canonical pulsars due to emission that is produced farther out in the magnetosphere. This is supported by recent studies with the NICER telescope, which show that MSPs radio profiles could originate in the outer edge of the beam instead of from the core of the emission beam (Guillot et al. 2019). Gentile et al. (2018) published fully-calibrated polarization profiles at 430 MHz, 1400 MHz, and 2300 MHz for 29 MSPs based on the NANOGrav 11-year data set (Arzoumanian et al. 2018) using the Arecibo Telescope. As expected, analysis of these profiles showed position angles that are generally inconsistent with the RVM. They also found microcomponents, which they defined as pulse components with peak intensities much lower than the total pulse peak intensity, in three pulsars.
In this paper, we present polarization profiles 1 of 23 MSPs observed with the Green Bank Telescope (GBT) at both 820 MHz and 1400 MHz as part of the NANOGrav 12.5-year data set (Alam et al. 2020). We 1 These data are available to be downloaded from data.nanograv.org/polarization. measure how the rotation and dispersion measures, and hence B (i.e. from Equation Eqn. 4), vary over the course of the data set.
In §2, we detail the observations. In §3 we discuss the polarimetric calibration, Faraday rotation fits, ionospheric corrections, and magnetic field calculations. In §4, we detail the results of the calibration and discuss the pulse profiles (comparison to published profiles, microcomponents, and frequency evolution/emission geometry), variations in dispersion measure and magnetic field, correlations with spindown parameters, polarization fractions, and implications for timing. In §5, we conclude the work.

OBSERVATIONS
We present a subset of the NANOGrav 12.5-year data set taken between MJDs 55265 and 56739 (2010 March 10 and 2014 March 23) at 820 MHz and 1500 MHz with the GUPPI instrument (DuPlain et al. 2008). We analyze observations of 23 pulsars, two of which overlap with Gentile et al. (2018) (PSRs J1713+0747 and B1937+21). Most pulsars were observed on a monthly cadence, with the exception of PSRs J1713+0747 and J1909−3744, which were observed weekly starting in 2013. The data were taken with the GUPPI backend at 820 MHz and 1.4 GHz with bandwidths of 200 MHz and 800 MHz, respectively. The data were coherently dedispersed, with frequency resolution of 1.56 MHz, and on average each observation lasted around 25 minutes. Table 1 shows the data timespan and number of observations for each pulsar at each frequency.
The data were run through the standard NANOGrav radio frequency interference (RFI) excision pipeline; for each frequency channel, the minimum and maximum values in the off-pulse region were found and any channels for which this value was an outlier relative to the surrounding channels were zapped (see Alam et al. 2020 for more details).
NANOGrav measures a DM at nearly every observation epoch. The value is then recorded as a DMX parameter in TEMPO, where DMX is the difference between the accepted DM and the measured DM at each epoch (Jones et al. 2017).
NANOGrav timing observations with the GUPPI data acquisition instrument began in 2010 March and continued into 2020. However, a technical problem arose in 2014 March, making all data collected after this date unsuitable for polarimetric work. The problem was instability in the time alignment of the digitizers for the X-and Y-polarizations of the telescope signal. This corrupted the polarization cross-products and made it impossible to recover full Stokes parameters from these Note-These numbers reflect only the data used in the analysis; the outliers have been removed.
data. The power in the two individual polarizations was uncorrupted, and well-calibrated total intensity measurements could still be derived, allowing for the use of these data in timing even without full Stokes parameter information. This instability only affects the polarization of the observations, should not affect the total intensity and therefore the timing after 2014.

Calibration Method
All NANOGrav observations go through a basic polarization calibration procedure. At the telescope, a 25-Hz broadband signal is generated at a noise diode and injected into the receiver. At the beginning of each observation, this artificial noise signal is split into two polarization signal paths and measured with the pulsar backend.
A calibration scan is taken for every NANOGrav observation. The noise signals themselves, and also the power in both X-and Y-polarizations, are calibrated by observations on and off a bright, unpolarized source (for the GBT, this is quasar B1442+101;NANOGrav Collaboration et al. 2015).
A set of four scans: pulsar, noise diode (which is the off-quasar scan), pulsar and noise diode, and quasar and noise diode, are used to obtain flux and polarization calibration solutions. A noise diode is observed with every pulsar scan but B1442+101 is observed once per each multi-day observing session at each frequency. This constitutes the standard calibration scheme, which is applied to all NANOGrav observations. While likely sufficient for timing purposes, in order to study the polarization in detail, more rigorous and precise polarization calibration is needed.
In this analysis, we used long-track observations of two pulsars to calculate Mueller matrix solutions. For our 820 MHz data, we used observations of PSR B1929+10 acquired by Kramer et al. (in prep.) for the double pul- sar, which were shared with NANOGrav. This pulsar is known for being very bright and has well-known polarization characteristics. We solved for the Mueller matrix at six epochs (MJDs 56244, 56419, 56608, 56793, 56984, and 57890) and used the solution closest to the epoch of each pulsar observation to calibrate the 820 MHz data. The solutions produced calibrated profiles for PSR B1929+10 that matched those in the literature (e.g., Stairs et al. 1999;Dai et al. 2015;Gentile et al. 2018) for every epoch, which suggests that our solutions accurately calibrated the data. See van Straten (2004) for full details of our calibration procedure.
At 1500 MHz, we used a single long-track observation of PSR J1022+1001 taken on MJD 55670 (2011 April 19) to calculate a Mueller matrix solution. Note that while PSR J1022+1001 has been found to show pulse profile variations by at most a few percent over the course of a year ), we do not expect this to affect our observations, as the solution was derived from an observation of this pulsar on a single day. After calibrating all of the data with the single solution, we found that the profiles were similar to both those in the literature (e.g., Stairs et al. 1999;Dai et al. 2015) and to each other, suggesting that using a single solution for multiple epochs produces accurately-calibrated profiles. Figure 1 shows an example of an 820 MHz solution used to calibrate our data. Panel (a) shows θ, the degree of cross-coupling between the receivers. Panel (b) shows ǫ, which indicates how much Stokes Q has leaked into Stokes V. The slight leakage of one Stokes parameter into another is caused by a small amount of nonorthogonality in the receivers. Panel (c) shows φ, the differential phase, which quantifies the mixing of the Stokes U and V parameters. Panel (d) shows γ, the differential gain. Ideally, γ = 0; in our data set, this parameter is consistent with zero for nearly all the epochs, with only slight offsets. Finally, Panel (e) shows G, the absolute gain for the receiver. As described earlier, we measured six independent realizations of the Mueller matrix as a function of frequency at 820 MHz at six different epochs. These realizations were generally consistent with each other.

Fitting for Faraday Rotation
To fit for Faraday rotation and calculate RMs, we used the rmfit feature of PSRCHIVE (van Straten et al. 2012), specifically the brute force method followed by the iterative position angle refinement technique. The iterative position angle refinement begins by using the brute force method to find the RM at which the linear polarization is maximized by fitting a Gaussian to the linear polarization vs. RM curve and using the centroid of the function as the best RM. We first re-binned the profiles to four frequency channels and 512 pulse phase bins at 820 MHz and 16 frequency channels and 512 pulse phase bins at 1500 MHz. We then searched in a range of -200 to 200 rad/m 2 with 200 steps for the majority of pulsars. An example of a fit is shown in Figure 2. Because of its location behind an HII region (Ocker et al. 2020), PSR J1643−1224 has a large RM ∼ −308 rad/m 2 (Yan et al. 2011a), so it requires finer frequency resolution to track the shift of PA with frequency. We therefore did not bin-or frequency-scrunch (which left us with the full 2048 bins and 512 channels) and searched from -550 to -150 rad/m 2 with 200 steps for the RM.
Once we had calculated an initial RM from the brute force method, we applied position angle refinement, which compares the position angles measured from the integrated profiles in the two halves of the band, and the weighted differential polarization angle (∆PA) is computed between the two halves of the band, using only the pulse phase bins in which the linear polarization is more than 3σ above the off-pulse noise. If ∆PA is larger than its uncertainty, the data are corrected for Faraday rotation using that RM and ∆PA between the two bands is estimated again. This process is repeated until ∆PA is smaller than its uncertainty, at which point the final RM is reported. This produces a more accurate RM estimate and uncertainty than the brute force method alone.
If an RM was not able to be fit with these parameters, we removed the profile from further analysis. For the most part, the number of observations taken out for this reason was relatively small (<15% of the total number of observations for each pulsar) but for PSRs J0645+5158, J0740+6620, J1455−3330, J1747−4036, and J1832−0836, the percentage removed was 25%, 43%, 32%, 31%, and 22%, respectively.
To ensure there were no outliers in RM values due to instrumental effects or miscalibration, we calculated the mean and RMS variations of the RMs for each pulsar and then removed data with RMs that were more than three standard deviations away from the mean from further analysis. After the first cut, a new mean was calculated and anything more than 3σ away from that value was cut. This process was repeated three times. Epochs with outlier RMs are not present in the combined (composite) profiles ( Figures 5-16) and were not used in the variability analyses. Most outliers showed up on specific days at both frequencies.
In addition to the method described above, we inspected the profiles by eye and eliminated any that looked noticeably different from the others. Criteria for this removal include incorrect handedness of the polarization, unusual variations in the profile baseline, and severe deviation from the composite profile on one epoch, all artifacts of a technical/instrumental problem with the observation.
Nearly all of the data sets that required outlier removal had <17% of observations removed. The exceptions were PSRs J0740+6620 (for which we excised 33% of the 1500 MHz data), J0931−1092 (which has 33% removed at 820 MHz), and J1832−0836 (which has 25% removed at 820 MHz). These high percentage are due to the small number of total observations relative to the number of excised observations. Though most outliers point to instrumental effects, high RMs that occur when a pulsar's line of sight passes close to the Sun may be due to a contribution from the solar magnetic field. We compared the epochs of the outliers we identified with those at which the relevant pulsar has the smallest elongation (the angle between the Sun and the pulsar). We also searched for outlier RMs at epochs at which DM peaks were detected. We find two such points for one pulsar, PSR J1614−2230, that are close to minimum elongation, when our line of sight to the pulsar passes closest to the Sun. See Section 4.2.2 for an in-depth discussion of these points.

Ionospheric Corrections
As the radio waves from the pulsar travel along our line of sight, they pass through the magnetic field of the Earth's ionosphere, which contributes a non-negligible amount to the measured RM. Therefore it must be subtracted in order to study the Galactic magnetic field. We used the ionFR (Sotomayor-Beltran et al. 2013) code, which uses publicly available GPS-derived total electron content CODE maps with a 2-hour time resolution for each day of observations, along with the eleventh release of the International Geomagnetic Reference Field, which covers the period when our data were taken. The code calculates the contribution of the ionosphere to the RM along the line of sight and takes into account the time of day of the observation, telescope location, and sky coordinates of the pulsar to get an accurate measurement for each hour of the day. We subtracted the ionospheric correction for the closest hour to the mid-point of each observation and were left with the RM due to the magnetic field of the ISM. Systematic uncertainties have been associated with this method, including a daily and yearly time dependence, with corrected RMs found to be accurate to 0.06-0.07 rad/m 2 (Porayko et al. 2019) . Therefore, we do not expect systematic uncertainties to be important for our RM measurements, given that the RM errors we derive are higher than this level (see Table 2).

Magnetic Field Calculations
To accurately constrain magnetic fields along the line of sight to the each pulsar (see Equation 4), we need to take into account variations in DM. The NANOGrav data set includes DMX parameters, which measure how much the DM of an observation varies from some fiducial or reference DM (see Jones et al. 2017). Table 2 shows the distance to each pulsar (calculated from parallax measurements from Alam et al. 2020), the reference DM (obtained from the par file for each pulsar), the average RM at each frequency (both corrected and uncorrected for the ionosphere), and the average magnetic field derived from the ionosphere-corrected RMs using Eqn. 4.
The uncertainties on the RMs show that the values are broadly consistent between the two frequencies, though some are discrepant at the 1 to 2-sigma level, suggesting that the error bars on the measurements are underestimated.
The error on the magnetic field at each epoch is the error from the RM and DM added in quadrature. The magnetic field value listed is the average over all epochs and both frequencies for each pulsar. Figure 3 shows the value of the magnetic field of pulsars around the sky using the values from this work combined with those of Gentile et al. (2018). The results are consistent with those of Sobey et al. (2019a), which uses pulsars and extragalactic sources in the northern sky to map the Faraday rotation measures, and hence the magnetic field of the Galaxy. For the most part, our results also match those of Gentile et al. (2018) as well as the values of Dike et al. (2020) which uses the Long Wavelength Array to analyze polarization of pulsars below 100 MHz.

Pulse Profiles
Here we present the polarization-calibrated average profiles from the method described in the previous section. Figures 5-16 show the composite profiles, which were made by summing the profiles from individual epochs. The position angle, which is shown in the top Table 2. Properties of each pulsar and derived quantities. The uncorrected RM at each frequency is the average of all of the measurements at each frequency, while the error is the standard deviation of all of the measurements divided by the square root of the number of measurements. The corrected RM at each frequency is the average of all of the measurements after the ionospheric correction is subtracted from each day, and the error is the standard deviation of all of those corrected measurements divided by the square root of the number of measurements. The magnetic field errors are a combination of the rmfit, ionFR, and DMX errors.
13.0 13.6 15.5 (7) 12.0 (6) 17 (1) 13 (1) 1.12 (6) (2) Note-The uncorrected RM at each frequency is the average of all of the measurements at each frequency, while the error is the standard deviation of all of the measurements divided by the square root of the number of measurements.
The corrected RM at each frequency is the average of all of the measurements after the ionospheric correction is subtracted from each day, and the error is the standard deviation of all of those corrected measurements divided by the square root of the number of measurements. The magnetic field errors are a combination of the rmfit, ionFR, and DMX errors. The quoted errors are the uncertainties on the last digit.
panel of each composite profile, is plotted when the linear polarization is >3σ above the off-pulse noise. Table  3 shows the fractions of total power of the emission, average fractional linear polarization, average fractional circular polarization, and average fractional absolute circular polarization of all pulsars in the data set at 820 MHz and 1500 MHz, all fractions are calculated with respect to the total power. Table 4 shows all previously published profiles for these pulsars. We present the first published polarization profiles at any frequency for PSRs J0636+5128, J0645+5158, J0740+6620, J0931−1902, J1125+7819, J1614−2230, J1747−4036, J1918−0642, and J2302+4442. We find no major discrepancies between our profiles and those previously published. The only exception is the sign of the circular polarization; Dai et al. (2015) uses the Institute of Electrical and Electronics Engineers definition of circular polarization whereas we use the IAU convention. This results in a sign change in the circular polarization (in the IEEE convention, left-hand circular polarization is positive and right-hand circular polarization is negative, whereas the IAU convention is the opposite).

Comparison to Published Polarization Profiles
Another exception is PSR B1937+21. At 820 MHz, the degree of linear polarization for B1937+21 shows epoch-to-epoch variability of up to ∼18%in the second main structure (the interpulse) . Because the RMs matched published values, we chose to carry out the analysis with them; the average profile is also similar to the literature, so we chose to present it. We will explore the reason for this variability in future work.
Overall, our RMs also agree with those previously published. There are several ways to measure RMs from pulsar profiles, and these methods have different systematic uncertainties. Most studies, such as this work and Yan et al. (2011b), use the rmfit method to calculate RMs and uncertainties, but other methods exist. For example, Sobey et al. (2019b) calculate Faraday RMs through Faraday spectra, or Faraday dispersion functions, with uncertainties calculated vis the method described in Brentjens & de Bruyn (2005). Our RM values are consistent within a few sigma of both of those results for the pulsars analyzed by both methods. The RM er-rors derived through these different methods are also consistent. This is reassuring, especially as Sobey et al. (2019b) calculate the RMs in a different way and using different bandwidths, center frequencies, and another telescope.

Microcomponents
We detect microcomponents in the pulse profiles of seven pulsars in this work. Microcomponents were discussed in Gentile et al. (2018) and here we define them as components that are <3% of the intensity of the highest peak on the average profile. Out of these seven, four pulsars have microcomponents that are detected for the first time. The microcomponents have varying degrees of polarization; for example, the microcomponents of PSR J2145−0750 are almost fully-polarized, whereas those of PSR J1909−3744 exhibit very little polarization. There is no apparent correlation between the amount of polarization in the microcomponents and that in the main pulse (i.e. the microcomponent of J2145−0740 is almost fully-polarized whereas the profile shows little).
Microcomponents that have been previously detected in other works have a flux density above 1.6 mJy, and all of the new ones have a flux density of less than 1.5 Table 3. Polarized intensity parameters. P is the phase-averaged power, L is the phase-averaged linear polarization, V is the phase-averaged circular polarization, V is the phase-averaged absolute value of the circular polarization, and I is the total intensity. The polarization fractions reported are those of the composite profiles. mJy. Because of our long data sets, which produce a very high S/N composite profiles, we are able to detect these very faint microcomponents. To ensure that the microcomponents were not an instrumental effect, we split each frequency band in half to see if the microcomponent was detected in each half. This was generally the case at both 820 MHz and 1500 MHz; the exception is J1713+0747, which exhibits a microcomponent only at 1500 MHz. This can be explained by the pulsar's very flat spectrum (Dai et al. 2015), resulting in lower S/N at lower frequencies.
The tests show that microcomponents are not an anomalous instrumental artifact but are of astrophysical origin. The detection of microcomponents demonstrates that MSPs emit over a wide phase range due to their larger opening angles and emission produced further out in the magnetosphere ).
These microcomponents make it difficult to define the duty cycle of millisecond pulsars. As noted in Gentile et al. (2018), they may cause an overestimation of the radiometer noise in the off-pulse region which could affect flux calibration (although NANOGrav does not rely on the radiometer noise for flux calibration). If these microcomponents are present in other pulsars, they would be revealed by longer data sets (and therefore higher S/N profiles). Microcomponents are generally most prevalent in our highest S/N pulsars; higher gain telescopes like the MeerKAT telescope in South Africa would improve that S/N, allowing us to probe weaker pulsars for these microcomponents (e.g., Spiewak et al. submitted). If not accounted for in template profiles, these mircocomponents could lead to higher uncertainties in TOA calculation. To make template profiles for TOAs, NANOGrav aligns and averages the reduced data profiles, and applies wavelet smoothing to the average profile (Alam et al. 2020). This wavelet smoothing preserves the microcomponents, and therefore they are taken into account when calculating TOAs.

Frequency Evolution/Emission Geometry
149 MHz (7) , 610 MHz ( The profiles for the majority of canonical pulsars are thought to evolve in frequency according to the core double cone model of Rankin (1983). This model makes specific predictions about how the number of components in a pulsar's average profile will vary with frequency. For example, a conal single pulsar will have two components at low frequencies (∼100 MHz) that will merge into one at higher frequencies (∼1 GHz). Xilouris et al. (1998) show that MSPs show three types of evolution: they can evolve minimally, evolve as predicted, or evolve contrary to any prediction. In their survey, 12 pulsars evolved minimally, five as predicted, and eight against predictions (e.g., with more components at higher frequency). This suggests that the emission of MSPs does not behave like the emission of canonical pulsars.
Frequency evolution is difficult to track in our pulsars, as many have more than five components and multiple structures in their profile (e.g. PSR J0931−1902). Out of the 22 MSPs for which we have accumulated profiles at both 820 MHz and 1500 MHz, 14 show the same number of components at both frequencies (i.e., develop min-imally) and eight seem to develop more components at higher frequencies, seemingly in contrast to the predictions of Rankin (1983) and in line with the Xilouris et al. (1998) results. While some of this evolution in MSPs could be due to decreased scatter broadening (causing separate components to appear as one at low frequencies), it supports the suggestion that MSPs do not evolve like canonical pulsars. While profiles evolve with frequency for all of the MSPs studied, there is no consistent trend and the frequency evolution is less dramatic than seen for non-recycled pulsars. Johnston et al. (2008) shows that in slow pulsars, the overall polarization fraction decreases as frequency increases, though some components can show an increase with frequency. This could be a consequence of a geometric process or involve orthogonal polarization modes. Overall, we find that the mean polarization fractions of linear and circular polarization do not show a clear trend with frequency. The exception is |V | /I; 12 MSPs in this study have higher |V | /I values at 1500 MHz, while only six have a higher fraction at 820 MHz and four have identical fractions at both frequencies. This shows a hint of a correlation, and this correlation is opposite to that than observed for canonical pulsars. However, this is a very small sample and further study is needed to confirm if this is the case in all millisecond pulsars.
As expected, many of the millisecond pulsars feature emission over a large portion of the profile (e.g. PSRs J1614-2230 and J2302+4442). This supports the idea that millisecond pulsar beams are wider than those of canonical pulsars due to emission produced farther out in the magnetosphere.
The position angle (PA) sweep is shown in the top panel of Figures 5-16; many of our pulsars (e.g. PSRs J1455−3330, J1918−0642, and J2145−0750), show very complex PA sweeps, which would require a model more sophisticated than the RVM. Only two pulsars in our data set show a quasi-S-shaped curve in the PA. Using PSRCHIVE, we searched an 18 by 18 grid in which α (the angle between the spin axis and the magnetic axis) and ζ (the sum of α and the angle between the magnetic axis and sightline β) are varied from 5 to 175 degrees in steps of 10 degrees. We perform this fitting for both pulsars at each frequency. The only significant result is for the L-band observation of PSR J1600−3053, where α = 162.8 ± 5.9, β = 2.35 ± 8.9 for a fit that has a χ 2 r value of 11.05. This shows that PAs are very difficult to fit in MSPs and a more sophisticated model incorporating emission far from the polar caps and/or more complex magnetic field structures is required to fit the position angle sweeps.

Variations in Measured Values
For each of the three parameters (ionosphere-corrected RM, DM, and B ), we performed a least-squares fit weighted by the uncertainties for a purely linear trend, a purely sinusoidal trend, and a combination of the two for all pulsars for which we have greater than one year of data. We only performed a sinusoidal and combination fit if a significant period with a false alarm probability (FAP) less than 5% was first identified through a Lomb-Scargle periodogram analysis. This FAP was calculated using the formula from Scargle (1982), which uses the length of the dataset and power spectral density to determine the probability that the period of the Lomb-Scargle periodogram is detected by random chance. That period was then used as the initial guess for the fitting. The reduced chi-squared (χ 2 r ) values were calculated for each fit; the trend reported for each pulsar is the model with the smallest χ 2 r value. The parameters for the trends are reported in Tables 4.1.3 and 6 and the data containing the best fit trend lines are shown in Figures 17−21. We plot the two fre-quencies separately in order to gauge which trends are truly astrophysical. In addition, one frequency may be more sensitive than another due to the pulsar's spectral index, DM, or RFI, so we may only see the trend significantly in one.
In the absence of astrophysical variations, we would expect the root-mean-square deviation of RMs to equal roughly the average 1-sigma error on those measurements. In Figure 4, following Caleb et al. (2019), we plot the ratio of the average RM error to the standard deviation vs. S/N. We find that these values are typically smaller than one, indicative of either real astrophysical variations or underestimated errors. We see more variation in RM values at higher S/N values. This seems to indicate that in the moderately high S/N regime RM errors are accurate, but that in the very high S/N regime, RM errors may be underestimated. Note that Caleb et al. (2019) found that RM errors measured for very low S/N profiles ( 17) were also underestimated. There are likely systematic effects that are not taken into account at high S/N, as shown by the lack of significant trends in the plots for bright pulsars such as J1713+0747.

DM Variations
Dispersion measure trends are shown in Table 4 There are twelve pulsars that overlap between our data sets, and Jones et al. (2017) finds significant trends in eleven, whereas we only find trends in three. For those pulsars in which we do both find significant trends, the slopes are roughly the same magnitude and the trends are the same for two of them (we find an extra sinusoidal trend in PSR J1012+5307). The differences in our results can be attributed to a lack of overlap in the data sets. Jones et al. (2017) is sensitive to longer term trends because they fit nine years of data, whereas we only include four years in our analysis. Yan et al. (2011a) also fit DM trends with one year of data observed at 1.4 GHz from You et al. (2007), though they fit only for linear trends. Their data are not sensitive enough for high-precision DM measurements, and they therefore only report upper limits on the slope of DM variations. For the pulsars that overlap between their paper and this one, the upper limits for only three Any period that had less than a 5% false alarm probability was not considered significant. The trend reported is the one with the smallest χ 2 r value. The pre-fit χ 2 r value refers to the χ 2 r of fitting a horizontal line through the data.  Table 6. Magnetic field trends. The results of fitting a linear trend, a purely sinusoidal, and a sinusoidal + linear trend to the magnetic fields. A weighted least-squares fitting routine was performed and the periods of the sinusoidal fits were first estimated with a Lomb-Scargle periodogram and then refined in the fitting routine. Any period that had less than a 5% false alarm probability was not considered significant. The trend reported is the one with the smallest χ 2 r value.The pre-fit χ 2 r value refers to the χ 2 r of fitting a horizontal line through the data. (PSRs J0613−0200, J1024−0719, and B1937+21) are significant. We find no significant trends in the DM of any of those pulsars. Discrepancies could be caused by our longer baselines. Though Yan et al. (2011a) predict that the slopes they measure are believed to be representative of the longer-term gradients, the linear trends they see are most likely fitted out over longer data sets (which is seen in our analysis). In addition, Donner et al. (2020) analyzed DM variations in 36 MSPs at a frequency of 150 MHz in a data set that spans 2012−2020. Nine pulsars overlap between our data sets. While they report linear trends for all nine of the pulsars, we find linear trends for only four of them. Their much lower observational frequency make them more sensitive to DM variations. This, combined with their longer datasets, likely explain this discrepancy; for the four pulsars for which we both measure trends, ours are generally of the same order of magnitude and are all of the same sign to those of Donner et al. (2020).

Variations in Measured B
The magnetic field variations are shown in Table 6 and the top panels of Figures 17-21. We find five pulsars with significant trends (J1600−3053, J1643−1224, J1713+0747, J1918−0642, and B1937+21). Four pulsars show a trend with a sinusoidal component, and three of the periods are consistent with one year, the other with a period of almost 700 days. Periods consistent with one year point to either contributions from the solar wind or magnetized clumps of material along our line of sight to the pulsar. We only see two to three full periods in the data set, so these are likely due to stochastic processes and are not true periodicities. As previously noted, because of corrupted data after the sampler board switch, we used a maximum of four years of data for each pulsar.
Two pulsars, PSRs J1643−1224 and J1918−0642, show significant linear trends in magnetic field. Assuming that this is due to movement along the line of sight through a region of increasing or decreasing Galactic magnetic field, we can use the timescale and slope  of the trends to calculate the ambient magnetic field over the distance the pulsars have traversed over the timespan of these observations. Local magnetic fields of roughly 400 mG for PSR J1643−1224 and 3200 mG for PSR J1918−0642 would be required over the distances of roughly 100 µpc traveled by the pulsars over the timespan of our observations in order to produce the changes in average magnetic field observed. This is much larger than ambient and/or local magnetic fields expected in the Milky Way.
van Ommen et al. (1997) measured the time variability of the RMs of PSRs B1556−44 and B1727−47 and found that local magnetic fields of 2 µG and 16 µG, respectively, were required to explain the observed variations. The latter was attributed to motion through irregularities within a nearby HII region. Rankin et al. (1988) observed RM and DM variations towards the Crab pulsar for two years and calculated a local magnetic field of ∼170 µG, consistent with its dense magnetic environment. Hamilton et al. (1985) observed another pulsar in a supernova remnant, the Vela pulsar, and found that the RM is increasing and found the magnetic field along the line of sight to be 22 µG, which was attributed to a magnetized cloud moving out of the line-of-sight to the source. Most recently, Johnston et al. (2021) used the ultra-wideband on the Parkes radio telescope to observe pulsars over two years. They measured the RM and DM and found that PSR J1825−1446 showed significant RM and DM changes, with the magnetic field along the line of sight changing by 0.2 µG in 2 years, which is due to the pulsar passing behind a magnetised filament in a supernova remnant.
Our ambient magnetic fields are much larger than any measured values, including the 1 mG fields sampled by PSR B1959−63 as it travels through the disk of its companion star (Johnston et al. 2005). This shows that the linear trends in magnetic fields we observe are much too large to be explained due to pulsar movement solely through an over-dense region along our line of sight, and are more likely due to our line of sight traversing variations in Galactic magnetic field structure in the transverse direction. Yan et al. (2011a) point out similarly large (∼0.1 mG) derived local magnetic fields for pulsars for which they measure linear changes in Galactic magnetic field with time (specifically PSRs J0613−0200, J1909−3744, and J2129−5721). Their slopes, however, are one to two orders of magnitude larger than ours. They use a different technique, relying on the slope of the RM divided by the slope of the DM to calculate the ambient magnetic field.
Their method, along with that of Hamilton et al. (1985), Rankin et al. (1988), andvan Ommen et al. (1997), assumes that the entire change in magnetic field is due to a small clump of material with a discrete RM and DM contribution into our line of sight, and does not account for the pulsar's movement along the line of sight. Our equation calculates the ambient magnetic field assuming that the magnetic field changes are due to the pulsar moving closer or further away from us in a region of dense magnetic field. If we make this as-sumption, the magnetic field for J1713+0747 (the only pulsar that shows a linear trend in both RM and DM) is 9 mG, which is more comparable to previous estimates but still large. Yan et al. (2011a) point to statistical fluctuations due to random spatial and temporal variations in the interstellar electron density and B to explain RM variations. Our numbers show that the magnetic field changes cannot entirely come from the motion of the pulsar through the interstellar medium. You et al. (2012) explored the effects of the Sun on pulsar RM values by observing PSR J1022+1001 when its line-of-sight passed close to the Sun. They found significant effects when the line-of-sight to the pulsar passed below 10 R ⊙ , which corresponds to ∼3 • of elongation.
We also checked the outliers for large changes in RM, DM, and B when the pulsars were close to minimum elongation. We found that PSR J1614−2230 experiences an increase in all three parameters when it came within 1.3 • of the Sun (which corresponds to ∼4.5 R ⊙ ). The increase in RM and DM at minimum elongation corresponds to a solar Galactic magnetic field contribution of 12(1) mG. This is consistent with You et al. (2012) and Ord et al. (2007), who report Galactic magnetic fields of the same order of magnitude at similar distances from the Sun.

Correlations with Pulsar Spin-Down Parameters
Studies such as Johnston & Kerr (2018) have examined the correlation between polarization fraction and spin-down parameters, but none have been conclusive.
Using the wealth of polarization information in this study, we examine the relationship between fractional linear and circular polarization and five parameters: spin period, age, surface dipole magnetic field, spin down energy loss, and proper motion. We find no conclusive evidence of any correlations between these parameters and linear or circular polarization fraction 820 MHz or 1500 MHz. If a relation did arise, it would give information about the magnetosphere, pointing to the fact that MSPs, for instance, with different spin periods have different sized magnetospheres. However, our sample size is fairly small, covering only a small range of distances, inclination angles, and other parameters. A larger sample size is needed for this analysis for any definitive conclusions to be drawn.

Timing Implications
Effects of polarization calibration on timing have been explored in many studies in the past decade, including Desvignes et al. (2016), Manchester et al. (2013), and Caballero et al. (2016). van Straten (2013) used matrix template matching to polarization calibrate PSR J1022+1001. They found that the RMS residuals decreased by a factor of two when polarization calibration was applied.
Pulsar time-of-arrival measurements calculated from data which have not been corrected for telescope polarization distortions, such as the NANOGrav 12.5-year data set (Alam et al. 2020), are susceptible to systematic timing uncertainties. These uncertainties will be higher for pulsars with larger polarization fractions (see Table 3). Correction of these data using the Mueller matrix formulation, as in the present paper, has the potential to improve the timing accuracy of such data sets. Also note that while incorrect polarization calibration could lead to higher levels of noise in the dataset, it would not show the spatial correlations expected for a gravitational wave signature. Rogers (2020) analyzed the effect of different combinations of polarization calibration meth-ods on the Parkes Pulsar Timing Array (PPTA) data using three techniques: Scalar Template Matching (STM), Measurement Equation Template Matching (METM), and Matrix Template Matching (MTM). STM, which is NANOGrav's method for calibrating profiles, models the transformation uses only the total intensity Stokes parameter. MTM, the method used in this paper, uses all four Stokes parameters to model the transformation between calibrated timing templates and the uncalibrated observations. METM, the method used by Gentile et al. (2018), relies on a bright pulsar as a standard source and produces a template/Mueller matrix solution for each day by forcing the observation of that pulsar to look like the template, obtaining a solution for each day. The work also relies on the Ideal Feed Assumption (IFA), which assumes that the receivers are perfectly orthogonal, the reference source is 100% polarized, and that the noise diode illuminates both receivers equally). Rogers (2020) calculated the TOAs for five millisecond pulsars using data calibrated with combinations of these techniques: IFA/STM, IFA/MTM, METM/MTM, and METM/STM. Both the METM combined with MTM and IFA combined with MTM method resulted in significantly more precise and accurate TOAs and timing residuals with smaller amounts of red and white noise, with the METM/MTM showing slightly better improvement overall. When compared to NANOGrav's method of IFA/STM, the combination of IFA/MTM used in Rogers (2020) improved the RMS of the post-fit residuals and the white noise residuals an average of 21% and 48% respectively, with a white noise residual improvement of above 60% in two pulsars.
Though METM produces TOAs with less red and white noise, it relies on the assumption that the pul-sars do not have any intrinsic polarization variability. It also removes any sensitivity to variability in the pulsars used as templates. However this work indicated the IFA/MTM method is just as effective as MTM/METM. Future work will apply the methods outlined in this paper to NANOGrav data to determine the effect of polarization-calibrated profiles on timing.

CONCLUSIONS
In this work, we presented polarization-calibrated profiles for 23 millisecond pulsars timed by the NANOGrav collaboration, which represent the first published polarization profiles for nine pulsars. NANOGrav's high S/N observations allowed for the discovery of very low intensity average profile components (microcomponents) in four pulsars. These are the highest S/N polarization profiles ever published for these millisecond pulsars and are made publicly available to the community to facilitate sensitive modeling of MSP emission mechanisms and geometries. We found that our MSPs are consistent with previous studies in that they evolve and behave differently than canonical pulsars.
We fit for Faraday rotation on each epoch and used the rotation measure and dispersion measure to calculate the magnetic field parallel to the line of sight of the pulsar. After fitting for a linear, sinusoidal, and sinudoisal + linear trend, we found a significant linear trend in three pulsars. Calculation of the ambient magnetic field produced large values on the order of micro-Gauss, which showed that the magnetic field changes cannot be entirely due to the motion of the pulsar along the line of sight and must be due to transverse motion through the large-scale Galactic magnetic field structure. Recent literature shows that this method of polarization calibration is likely to greatly improve the timing precision of our pulsars, which will be examined in future work.
These data only represent a portion of those obtained by the NANOGrav timing campaign. New ultra-wideband receivers on the GBT will provide more sensitivity. Also, the Canadian HI Mapping Experiment (CHIME) telescope will provide complementary frequency coverage to track how the polarization and microcomponents behave at lower frequencies. Author Contributions: HMW carried out the analysis and prepared the text, figures, and tables. MAM helped with the development of the framework and the text. PAG helped with the development of the framework. MLJ helped with the analyses and code. RS assisted with the analysis and the discussion. All authors contributed to the collection and analysis of the NANOGrav 12.5-year data set; see Alam et al. (2020) for further details.
Software: numpy (Oliphant 2006;Van Der Walt et al. 2011), matplotlib (Hunter 2007), PSRCHIVE  You, X. P., Coles, W. A., Hobbs, G. B., & Manchester, R. N. 2012, MNRAS, 422, 1160, doi: 10.1111/j.1365-2966.2012.20688.x You, X. P., Hobbs, G., Coles, W. A., et al. 2007, MNRAS, 378, 493, doi: 10.1111/j.1365-2966.2007.11617.x Figure 5. Pulse profile for pulsars J0340+4130 and J0613-0200. The black line is the total intensity, red is the circular polarization, and blue is the circular polarization. The polarization position angle is shown in the top panel. Figure 6. The pulse profile for pulsars J0636+5128 and J0645+5158. The black line is the total intensity, red is the linear polarization, and blue is the circular polarization. The polarization position angle is shown in the top panel. Figure 7. Pulse profiles for pulsars J0740+6620 and J0931-1902. The black line is the total intensity, red is the linear polarization, and blue is the circular polarization. The polarization position angle is shown in the top panel. Figure 8. Pulse profiles for pulsars J1012+5307 and J1024-0719 including microcomponents. The black line is the total intensity, red is the linear polarization, and blue is the circular polarization. The black arrow points to the location of the microcomponent in J1024-0719. The polarization position angle is shown in the top panel. The microcomponent plots for J1024-0719 have been plotted with fewer bins to increase the signal-to-noise. The polarization position angle is shown in the top panel. Figure 9. Pulse profiles for pulsars J1125+7819 and J1455-3330 including microcomponents. The black line is the total intensity, red is the linear polarization, and blue is the circular polarization. The black arrow points to the location of the microcomponent of J1455-3330. The microcomponent plots for J1455-3330 have been plotted with fewer bins increase the signal-to-noise. The polarization position angle is shown in the top panel. Figure 10. Pulse profiles for pulsars J1600-3053 and J1614-3053 including microcomponents. The black line is the total intensity, red is the linear polarization, and blue is the circular polarization. The black arrow points to the location of the microcomponent of J1600-3053. The microcomponent plots for J1455-3330 have been plotted with fewer bins increase the signal-to-noise. The polarization position angle is shown in the top panel. Figure 11. Pulse profiles for pulsars J1643-1224 and J1713+0747 including microcomponents. The black line is the total intensity, red is the linear polarization, and blue is the circular polarization. The black arrow points to the location of the microcomponent in each J1713+0747 profile. The microcomponent plot for J1713-0747 have been plotted with fewer bins to increase the signal-to-noise. The polarization position angle is shown in the top panel. Note: there is no detection of the microcomponent of J1713+0747 at 820 MHz, the plot is just shown for comparison. Figure 12. Pulse profiles for pulsars J1744-1134 and J1747-4036. The black line is the total intensity, red is the linear polarization, and blue is the circular polarization. The polarization position angle is shown in the top panel. Figure 13. Pulse profiles for pulsars J1832-0836 and J1909-3744 including microcomponents. The black line is the total intensity, red is the linear polarization, and blue is the circular polarization. The black arrow points to the location of the microcomponent in each J1909-3744 profile. The polarization position angle is shown in the top panel. Figure 14. Pulse profile for pulsars J1918-0642 and B1937+21 including microcomponents. The black arrow points to the location of the microcomponent in each B1937+21 profile. The black line is the total intensity, red is the linear polarization, and blue is the circular polarization. The polarization position angle is shown in the top panel. Figure 15. Pulse profiles for pulsars J2010-1323 and J2145-0750 including microcomponents. The black arrow points to the location of the microcomponent in each J2145-0750 profile. The black line is the total intensity, red is the linear polarization, and blue is the circular polarization. The polarization position angle is shown in the top panel. Figure 16. Pulse profiles for pulsars J2302+4442. The black line is the total intensity, red is the linear polarization, and blue is the circular polarization. The polarization position angle is shown in the top panel. Figure 17. Dispersion measure and magnetic field changes over time for pulsars J0340+4130, J0613-0200, J0645+5158, and J1012+5307. The uncertainties on the DM come from those on the DMX value and uncertainties on the magnetic field are a combination of the uncertainties on the ionosphere-corrected RM (which are a combination of those of fitting for Faraday rotation and from the ionospheric correction) and the DM. No trendlines are shown because the lowest χ 2 r value for the fits was that of a horizontal line with a slope of zero.  Figure 18. Dispersion measure and magnetic field variations over time for pulsars J1024−0719, J1455−3330, J1600−3053, and J1614−2230. The uncertainties on the DM come from those on the DMX value and uncertainties on the magnetic field are a combination of the uncertainties on the ionosphere-corrected RM (which are a combination of those of fitting for Faraday rotation and from the ionospheric correction) and the DM. Any trendlines shown represent the trend with the lowest χ 2 r value. If no trendlines are shown then the lowest χ 2 r value for the fits was that of a horizontal line with a slope of zero. Note: the plots for J1614-2203 contain two outliers at epochs of small ecliptic angle (less than 3 degrees) (MJDs 55892 and 55893, as discussed in Section 4.2.2). These points are excluded from the fitting and the mean RM and B calculation but included in the plot to show the spike in RM, DM, and B when the pulsar is close to the Sun.  Figure 19. Dispersion measure and magnetic field changes over time for pulsars J1643−1224, J1713+0747, J1744−1134, and J1747−4036.The uncertainties on the DM come from those on the DMX value and uncertainties on the magnetic field are a combination of the uncertainties on the ionosphere-corrected RM (which are a combination of those of fitting for Faraday rotation and from the ionospheric correction) and the DM. Any trendlines shown represent the trend with the lowest χ 2 r value. If no trendlines are shown then the lowest χ 2 r value for the fits was that of a horizontal line with a slope of zero.  Figure 20. Ionosphere-corrected rotation measure, dispersion measure, and magnetic field changes over time for pulsars J1909−3744, J1918−0642, B1937+21, and J2010−1323. The uncertainties on the DM come from those on the DMX value and uncertainties on the magnetic field are a combination of the uncertainties on the ionosphere-corrected RM (which are a combination of those of fitting for Faraday rotation and from the ionospheric correction) and the DM. Any trendlines shown represent the trend with the lowest χ 2 r value. If no trendlines are shown then the lowest χ 2 r value for the fits was that of a horizontal line with a slope of zero.  Figure 21. Dispersion measure, and magnetic field changes over time for pulsars J2145−0750 and J2302+4442. The uncertainties on the DM come from those on the DMX value and uncertainties on the magnetic field are a combination of the uncertainties on the ionosphere-corrected RM (which are a combination of those of fitting for Faraday rotation and from the ionospheric correction) and the DM. No trendlines are shown because the lowest χ 2 r value for the fits was that of a horizontal line with a slope of zero.