The Thousand-Pulsar-Array programme on MeerKAT XIII: Timing, flux density, rotation measure and dispersion measure timeseries of 597 pulsars

We report here on the timing of 597 pulsars over the last four years with the MeerKAT telescope. We provide Times-of-Arrival, pulsar ephemeris files and per-epoch measurements of the flux density, dispersion measure (DM) and rotation measure (RM) for each pulsar. In addition we use a Gaussian process to model the timing residuals to measure the spin frequency derivative at each epoch. We also report the detection of 11 glitches in 9 individual pulsars. We find significant DM and RM variations in 87 and 76 pulsars respectively. We find that the DM variations scale approximately linearly with DM, which is broadly in agreement with models of the ionised interstellar medium. The observed RM variations seem largely independent of DM, which may suggest that the RM variations are dominated by variations in the interstellar magnetic field on the line of sight, rather than varying electron density. We also find that normal pulsars have around 5 times greater amplitude of DM variability compared to millisecond pulsars, and surmise that this is due to the known difference in their velocity distributions.


INTRODUCTION
The radio emission from pulsars is thought to originate from close to the surface of the star and locked to its rotation.The radio emission is usually highly linearly polarised and is beamed along the magnetic axis, resulting in pulses of emission being detected by a distant observer.The time-of-arrival of the pulses at an Earth-based observatory yields much information, it tells us the rotation rate of the pulsar, how fast it slows down, how irregular the slow-down is and about transient events such as glitches.Together, these inform us about the evolution of pulsars over time and provides insight into the composition and structure of the neutron star interior.Additionally the traverse of the emission through the interstellar medium (ISM) gives information on the composition, structure and magnetic fields of the ISM via measurements of dispersion measure (DM, e.g.Yao et al. 2017) and rotation measure (RM, e.g.Han et al. 2018).
The Thousand Pulsar Array (TPA) programme on the MeerKAT telescope (Johnston et al. 2020) forms part of the larger MeerTime key science project (Bailes et al. 2020).The top-level goals of the TPA are two-fold.First to use the full sensitivity of the MeerKAT telescope ★ E-mail: mkeith@pulsarastronomy.net to provide a census of both the integrated and single-pulse properties for more than 1000 pulsars (Posselt et al. 2021;Johnston et al. 2023;Posselt et al. 2023;Song et al. 2023).The second major goal is to investigate the time-varying properties of pulsars through regular, high fidelity (Song et al. 2021) observations of a large (∼ 600) and diverse sample of pulsars.These observations have a median cadence of 27 days, though the exact observing pattern has changed several times over the course of the programme primarily due to scheduling constraints on the telescope.This paper describes the data products from the TPA regular monitoring programme that are published alongside this paper.This first data release consists of time of arrivals (ToAs) and pulsar ephemerdies for 597 pulsars, as well as measurements of spin-down rate, flux density, rotation measure and dispersion measure for each pulsar over the period 2019 March to 2023 May.Pulse profile shapes are also monitored within the TPA, details of which will appear in a subsequent paper, with initial results published in Basu et al. (2023;submitted).Section 2 briefly outlines the observations and defines the sample of pulsars in the data release.Section 3 describes the methods used to measure each of the parameters presented in the tables and figures comprising the data release, which are themselves described in Section 4. Some broad highlights of the data are presented in Section 5, with a particular focus on the measurements of DM and RM variations as an example of what may be done using the data.

Observations
Observations for the TPA programme fall into two categories.For the census projects, observations were carried out with at least 60 antennas of the MeerKAT array.Observations were done with the Lband receiver and, generally, every pulsar was observed once with an integration time determined using the prescription outlined in Song et al. (2021).For the monitoring project we aim to observe some 500 pulsars once per month.These observations split the array into two sub-arrays each consisting of ∼ 30 antennas.Here, observation times are very short (typically less than 2 mins in duration) in order to observe as many pulsars as possible in the monthly allocation of 16 hours of telescope time.The data in this data release combine both types of observations, though the vast majority of data values come from the regular observations.Although we aim for continuous observations with a regular sampling, around half of the pulsars exhibit a gap in observations due to periods where pressure on available telescope time required us to reduce the number of pulsars observed.Hence, there is a gap of around one year, from 2021 March to 2022 February, for 294 pulsars, and a gap of around 200 days, from 2019 Nov to 2021 May, for 74 pulsars.
A description of the observing system can be found in Bailes et al. (2020) and Johnston et al. (2020).In brief, we used the observational band from 896 to 1671 MHz with 928 frequency channels.The data are folded into sub-integrations each of length 8 s for the duration of the observation and there are 1024 phase bins per pulse period.Polarization calibration follows the method outlined in Serylak et al. (2021) and flux calibration is achieved as described in Posselt et al. (2023).Data outputs are in psrfits format (Hotan et al. 2004).In total, more than 24,000 observations have been made for the TPA programme over the 4-yr timespan.

The pulsar sample
There are in excess of 1250 pulsars observed as part of the TPA (Posselt et al. 2023;Song et al. 2023).Many of these only have a one or two census-style observations, but there are 614 pulsars that had at least 10 observations in the period from 2019 March to 2023 May.Pulsars were selected for the regular monitoring campaign based on the quality of data that could be obtained in the very short observation time, and the list was refined several times with the priority on ensuring data of good quality whilst also trying to keep a good spread of pulsars over the - diagram (see Figure 1) The latter criteria was included to reduce selection bias when the data are used for studies of properties of the pulsar population, though we note that the TPA does not observe millisecond pulsars (MSPs) and the MeerTime programme to observe MSPs is described elsewhere (Spiewak et al. 2022).From these 614 we exclude the 7 pulsars which have binary companions (PSRs J0045−7319, J0823+0159, J1141−6545, J1302−6350, J1740−3052, J1906+0746 and J1930−1852) as these require more careful analysis and will be described elsewhere.We further exclude three pulsars (PSRs J1717−4054, J1906+0414 and J1929+1357) which have a high nulling fraction and hence very few detections, and one pulsar (PSR J1644−4559) which was used mainly for system tests.Finally, there are six pulsars (PSRs J0835−4510, J1301−6305, J1420−6048, J1614−5048, J1718−3825 and J1907+0631) which we exclude as we were unable to obtain a coherent timing solution because of the presence of large glitches.There are therefore a total of 597 pulsars presented in this paper.

Templates and ToA measurement
ToA measurements are one of the data products from the the Meer-Pipe pipeline which runs automatically on all TPA observations.The ToAs are computed on the 8 sub-band, fully time-averaged, data using psrchive to cross correlate the observed pulse profile with a noise-free template, using the Fourier-domain markov-chain (FDM) method in pat (Hotan et al. 2004).A single template is used across all 8 sub-bands, generated from the 'best' TPA observation, which is typically the observation from the TPA census, which uses the full MeerKAT array and much longer intergration time than the typical TPA timing observations, or from averaging several observations where signal-to-noise ratio is low.For most pulsars the noise-free template is generated by means of a Gaussian process (GP) modelling of the observed average pulse profile (Johnston & Karastergiou 2019).To improve the performance of the template, the profile is tapered off to zero outside the on-pulse window computed in Posselt et al. (2021), which ensures that no information from the off pulse is used in the template matching.Each automated template was inspected and for 58 pulsars, notably those with scattering tails, interpulses or very complex shapes, the automated template generation was deemed unsuitable.For these pulsars, templates were computed by fitting von Mises functions1 'by hand' using paas from psrchive, or in the case of two highly scattered pulsars (PSRs J1701−4533 and J1828−1101), using von Mises functions convolved with a single exponential with a bespoke routine.Since only a single template was used across the whole band, users should be aware that the resulting ToAs will include effects of profile evolution which must be accounted for in any later timing analysis.
Initial timing was performed using tempo2 to compute and fit the timing residuals.The initial ephemeris was taken from psrcat, and coarsely refined whilst searching for possible phase wraps between observing epochs.For some pulsars with significant spin noise or glitch activity, the catalogue ephemerides lost phase at a rate of up to several turns per day, so for some cases we relied on ephemerides provided by long-term timing programmes with the Lovell telescope at Jodrell Bank Observatory, or the Murriyang telescope at Parkes Observatory.
The MeerPipe pipeline produces ToAs from all observing epochs, even if the pulsar was not detected, or the observation was corrupted by interference.Therefore some removal of bad data is needed, and ToAs were discarded if the pulsar was not detected, or if the ToA uncertainty was so great as to not provide useful information to the timing.In most cases the lack of detection is due to pulsars with known nulling, or weak modes that are not detectable in the typical TPA timing observations with less than 2 minute durations.

The timing model
The final TPA timing data product was refined from the coarse timing model by means of applying a Bayesian noise modelling simultaneously with the timing model fitting.The pulse noise model consists of a time-correlated achromatic spin noise, plus a time uncorrelated white noise component.The achromatic spin noise is implemented using the Fourier-domain GP as described in Lentati et al. (2013) to model the noise by a power-law parameterised by dimensionless amplitude  red and slope , and power spectral density given by where  yr is a frequency of 1 per year.The white noise consists of the widely used EFAC and EQUAD parameters that scale and add in quadrature the independent white noise for each observation, plus an additional ECORR term correlated across observations at each epoch.
For ToAs with formal error given by   , the white noise covariance matrix is given by where  is a diagonal matrix with Σ  =  2  ,  is the identity matrix, and  is a block-diagonal matrix such that Δ ,  = 1 for ,  from the same observing epoch and 0 otherwise (see e.g.van Haasteren & Vallisneri 2014 for a detailed description).
In addition to fitting for the noise parameters, we also marginalise over the following parameters in the pulsar ephemeris using the linearised model of tempo2, • F0, F1 and F2 -Spin frequency and two derivatives.
• RAJ and DECJ -the position of the pulsar in J2000 coordinates.
• Arbitrary phase jumps for 6 of 8 sub-bands.
The arbitrary phase jumps between sub-bands absorb any evolution of the pulse profile with frequency.In principle we may wish to fit for phase jumps between all 8 sub-bands, however, this would be fully covariant with the 'zero phase' parameter that is always fit within tempo2.Additionally, although the change in DM from epoch to epoch would be well measured, the average of the DM parameters would be unconstrained and hence not produce meaningful results.Therefore, we must leave two of the sub-bands phase without a phase jump in order to define the long-term average DM, and so that the fit is not degenerate with the zero phase parameter.In practice then the 'zero phase' parameter becomes defined by the average phase of these two sub-bands, and the long-term average DM is defined by the phase difference between these two sub-bands.In this analysis, we chose to keep fixed the top and bottom sub-bands, so that the average DM is determined using the most widely separated frequency channels, i.e. those centred at approximately 944 and 1623 MHz.We note that since the phase jumps are constant in time, the choice of sub-bands to fit only affects the average DM, and not the epoch-to-epoch variability, which is the main focus of this work.
The noise modelling is performed using run_enterprise (Keith & Niţu 2023;Keith et al. 2022), which is based on the enterprise framework (Ellis et al. 2019), and sampling is performed using emcee (Foreman-Mackey et al. 2013).The resulting noise model is then used to fit for the maximum likelihood values for the pulsar timing parameters using tempo2.

Timeseries of the spin-down rate
The timing residuals obtained for each pulsar are used to drive a time-domain GP, following Section 3.2 of Brook et al. (2016).We use the same software and process used for the recent analysis of PSR J0738−4042 by Lower et al. (2023) to obtain  and its associated uncertainty.The GP is constructed with a single squared exponential kernel and a white noise term, and the second derivative of the GP can be directly computed to derive ().We sampled the GP hyperparameter posterior distributions using bilby (Ashton et al. 2019) as a wrapper for the dynesty nested sampling algorithm (Speagle 2020).We then generated  points at each observing epoch directly from the GP.As always with this analysis, the assumption is that the residuals can be entirely explained by a time-variable .In the majority of cases, the GP models the residuals extremely well, but there may be examples where the use of an additional kernel, together with strong priors on the hyperparameters, would result in a better fit.We also caution that the uncertainties for  can depend on the choice of kernel, and are only correct under the above assumptions.

Dispersion Measure versus time
The DM of a pulsar measures the integral of the electron density,   along the line of sight, , (3) Observationally, the DM can be obtained by measuring the time delay between the arrival times ( 1 and  2 ) of the radio pulse at two different frequencies ( 1 and  2 ) via where  is the dispersion constant and is equal to 2.410 × 10 −4 MHz −2 cm −3 pc s −1 (see Kulkarni 2020 for a discussion of this value).In practice the DM to a pulsar is not a constant, and we measure the DM for each epoch as part of the pulsar timing process described in Section 3.1.2.It is worth noting that although the time dependence of the DM is usually well measured, there may be a systematic error in the DM if there is significant pulse profile evolution across the band, including effects from scattering.In order to get an indication of the scale of this systematic error, we compute the change in the estimated mean DM if we reference the DM against the highest or lowest pair of sub-bands rather than the most widely spaced sub-bands.It should be noted that this systematic error shifts all DM values up or down together, and does not change the perceived time variability.

Rotation Measure versus time
The rotation measure (RM) of a pulsar measures the integral of the product of the electron density and the magnetic field parallel to the line of sight where   and  are the mass and charge of an electron respectively.Observationally, the RM can be obtained by measuring the change in the position angle (PA) of the linearly polarized radiation at two observing wavelengths ( 1 and  2 ).
In practice, several different methods are employed to determine the RM and its error bar.
Here, the methodology as described in Ilie et al. (2019), and implemented in psrsalsa 2 (Weltevrede 2016) is used.It is based on the RM synthesis technique (RMST; e.g.Brentjens & de Bruyn 2005), which optimises the degree of linear polarization as function of RM.Uncertainties are estimated via bootstrapping the data.This relies on the identification of an on-and off-pulse region, which is done automatically using the same methodology as explained in Song et al. (2023), and relies on the profile templates (see also Section 3.1.1).
The RM as measured at the telescope is determined for each epoch, resulting in a timeseries of RM measurements.Observations with very weak detections of the pulse profile, for instance because of RFI, were ignored in the analysis.This was done by visually inspecting pulse profiles of the in total 24,592 observations.To make this process more efficient, the pulse profiles of each pulsar are rankordered by the signal to noise ratio (S/N).The pulse profiles for each pulsar were visually inspected starting with the lowest S/N, and observations were rejected until a pulse profile of sufficient quality was encountered.Higher S/N detections were accepted without further visual inspection.In total 703 observations out of a total of 24,592 are excluded from the RM timeseries in this way.
The Earth's ionosphere contributes to the value of RM and this varies with time.This contribution was predicted using the ionFR package3 (Sotomayor-Beltran et al. 2013), which relies on global Total Electron Content (TEC) values4 .The predicted ionospheric contribution to the RM is subtracted in the figures (such as Fig. 7, and the figures in the online material).In the tables of the online material the measured RM (which includes a contribution of the ionosphere) and the predicted RM contribution by the ionosphere are reported separately for each observation.

Flux density
Flux calibration is carried out following the procedure described in Posselt et al. (2023).In brief, the flux densities are scaled from the system counts via the radiometer equation and by using observations of pulsars at high galactic latitude (where the sky temperature is well understood) to calibrate the the sky temperatures of all other pulsars.The scatter in the observed flux density is generally much larger than the formal uncertainty.Intrinsically, the flux density of a pulsar is observed to be largely stable with a low modulation index (see e.g.Kumamoto et al. 2021).However, this is only true for observations of long duration (several thousand pulse rotations).In the monitoring program, a pulsar is typically observed for only 90 s or a few hundred rotations at best.The measured flux density will then fluctuate about some mean, depending on the random selection of single pulse flux densities from the underlying (and unknown) distribution.Furthermore, the traverse of the radio emission through the turbulent interstellar medium causes diffractive and refractive scintillation, which manifests itself as frequency and time dependent flux density variations.If these variations have comparable bandwidth and timescales to the observations, high modulation will be seen in the light curves (Kumamoto et al. 2021).Finally, pulsars are observed to have a steep spectral index (Jankowski et al. 2018;Swainston et al. 2022;Posselt et al. 2023) and averaging in frequency over a large bandwidth biases the flux densities towards lower frequencies than the nominal centre of the band.With all these caveats in mind, here we simply average the pulse profile data in time and frequency for a given observation and present a single flux density for that observation.

DESCRIPTION OF TABLES AND FIGURES
The data described in this section can be obtained from Zenodo, doi.org/10.5281/zenodo.8430591.

Timing files
Pulsar rotational ephemerides and ToAs are provided in tempo2 compatible 'par' and 'tim' files.The ToAs in the tim files are annotated with a number of flags from the meerpipe reduction pipeline, including signal-to-noise ratio (-snr) and observing time span (-length).Each ToA is also associated with a relative pulse number to enable phase tracking in tempo2.Each pulsar is supplied with a pair of par files.The first par file contains the complete output of the noise modelling, as produced by the Bayesian pipeline.We also provide a simplified par file without the noise modelling and DM measurements specific to our data that is more easily transferred to other software and/or combined with data from other telescopes.

Observation tables
We provide a table containing the list of 597 pulsars we are monitoring, along with the MJD of the first and last observation and the number of data points.This table also includes the mean DM as well as the systematic offset when using either the two highest bands or the two lowest bands to define the DM.For each pulsar, in addition to the par and tim files described above, we also supply 6 ascii tables, one each for the residuals, dispersion measure, rotation measure, the rotation frequency derivative ( ) and flux density versus observation epoch.Each file has a one-line header with column labels that are described in Table 1.Measurements that are obviously spurious have been removed where, for instance, visual inspection of the observed profile suggested that the pulsar was not detected with a sufficient signal-to-noise ratio to make a measurement of a particular parameter.Hence the number of rows in each table may differ due to the different sensitivity to each measured parameter.

Figures
For each pulsar we provide a visual representation of the data described above.An example plot is given in Figure 2 for PSR J1117−6154.The five panels show the timing residual (see

The interstellar medium
The strength of turbulence in the ISM can be characterised by a power-law over many orders of magnitude bounded by the inner and outer scales (Rickett et al. 1984).The various scales can be explored by scintillation observations (both diffractive and refractive), DM variations with time and RM variations with time.Given dense samples and a long time-baseline, the structure function is the implement of choice for exploring the nature of the turbulence (e.g.Cordes & Rickett 1998).However, the structure function is the average power as a function of lag, and requires being able to derive multiple independent samples of the process (at the lag of interest) from the data.For most pulsars in our sample, we are only sensitive to power on the timescale of our dataset, and hence the structure functions contains little meaningful data.We therefore resort to simply measuring any gradient in the DM versus time data, which gives an indication of the strength of the DM variations on the longest timescales.Table A1 records those pulsars for which we measure a gradient greater than 3- from zero.There are 87 pulsars with a significant gradient (out of 597), a significantly larger fraction than the 4 out of 160 pulsars reported in Petroff et al. (2013), highlighting the improved sensitivity in the TPA data.
One also expects that the RM should change with time, particularly if the DM is also changing.In addition, however, changes in the average magnetic field strength along the line of sight should also be observed, particularly if the pulsar is in or behind a magnetic structure such as a supernova remnant.In Table A2 we list those pulsars for which the slope of RM versus time is more than 3- significant.There are 58 such pulsars of which 15 are in common with pulsars with a DM gradient.
There are two pulsars which show very large changes in both DM and RM as shown in Figure 3. PSR J1114−6100 has a very large RM for its location and was discussed extensively in Johnston et al. (2021).Over the course of the monitoring campaign its RM has changed by some 6 rad m −2 and its DM by 1 cm −3 pc.PSR J1833−0827 had a large negative DM slope in Petroff et al. (2013).Here we see that the DM slope is positive and that the RM has changed by 15 rad m −2 in less than 2 yr.As Petroff et al. (2013) pointed out, this pulsar has an associated X-ray pulsar wind nebula (Esposito et al. 2011), making the nebula the most likely source of the fluctuations in RM and DM.

DM slopes vs DM
If the DM slopes we observe in our sample originate in discrete structures in the ionised ISM (IISM) then we may expect that the amplitude would scale with the square-root of the number of such structures encountered, and hence roughly with the square-root of the DM.Alternatively, we may expect that the observed DM slopes come from power on the longest length scales in the turbulent IISM, which are expected to be at least 100 pc (Armstrong et al. 1981), far longer than the scales probed by our ∼ 4 year dataset.To get a feeling for how the DM variations due to IISM turbulence may be expected to scale with DM, we can look at scattering, which is another effect of IISM turbulence.The 'level of turbulence' in the IISM can be estimated from measurements of interstellar scattering (Cordes 1986), where the scattering measure (SM) is the path-integrated variance in the electron density along the line of sight (Backer et al. 1993).The SM can be estimated from observations of scattering timescale,  scatt , where for Kolmogorov turbulence, SM ∝ ( scatt /) 5/6 , ( where  is the distance to the pulsar (Cordes et al. 1991).If the power in the IISM is independent of location then one expects to find  scatt ∝ DM 2.2 , however measurements of  scatt show evidence for inhomogenities in the level of turbulence, with Krishnakumar et al. (2015) finding over a the full range of DM values observed in the pulsar population.
Whilst the SM is a measure of the variance of the electron density, the amplitude of the DM variations should scale with the standard deviation of the DM, and hence all else being equal, the typical amplitude of DM slope to a pulsar should scale with the square-root of the expected SM (Backer et al. 1993).Hence, we can insert Equation 8 into Equation 7 and square-root.Making also the simplifying assumption that  ∝ DM we can get a scaling relation Hence we may expect that the exponent of DM slope to evolve slightly with DM, from 1/2 at low DM, to 4/3 at high DM, but in practice the intrinsic scatter means that this does not seem distinguishable from a DM slope linearly proportional to DM over the TPA sample.
Indeed, we find that the measured DM gradient for the TPA pulsars appears consistent with linear increase with DM with a slope of the order of 10 −4 ×DM per year.Figure 4 shows this for the TPA pulsars, and also for the sample of MSPs from the MeerTime PTA data release (Miles et al. 2023) which is much more sensitive to DM variations due to the nature of precision MSP timing.Although both seem to scale with DM, the TPA sample shows a factor of around 5 times larger DM slopes for a given DM.PTA pulsars are, due to selection effects in the discovery of pulsars, more likely to be at lower DM than the TPA sample, and hence some of the observed difference may be due to the inhomogenties in the IISM such as those which give rise to Equation 9.However, as seen in the dot-dashed line in Figure 4, the increased DM variation at large DM does not seem strong enough to explain the difference.
An important consideration for interpreting these results is that the rate at which we observe DM fluctuations does not only depend on the scale of the turbulence in the IISM, but also on the rate at which the line of sight to the pulsar traverses the IISM.For most pulsars, the line of sight velocity is dominated by the velocity of the pulsar, therefore, we may suppose that the the difference in scale of DM gradients may instead be due to the smaller velocity distribution observed in MSPs compared to the general pulsar population (Hobbs et al. 2005;Shamohammadi et al. 2024).Figure 5 shows the fractional DM slope as a function of an estimate of the pulsar 2-D velocity.For the PTA pulsar sample, the velocities are taken from (Shamohammadi et al. 2024), which are based on VLBI or timing parallax and proper motions.This velocity for TPA pulsars is estimated from proper motion and distance recorded in psrcat.Parallax measurements are used if available, or otherwise distance is inferred from the DM using the 'YMW16' electron density model (Yao et al. 2017).It should be noted that the uncertainties in the velocity of TPA pulsars are not rigorous as distance estimates can be unreliable, and we have made an assumption of normally distributed independent uncertainties on all measurements.We refer readers to Verbunt et al. (2017) or Shamohammadi et al. (2024) for a complete and careful study of the intrinsic distribution of pulsar velocities.Nevertheless, studies of pulsar velocities seem clear that the 'normal' pulsar population sampled by the TPA has an average velocity at least 3 times greater than the PTA sample, and there is evidence that a large fraction of the normal pulsars have velocities 5-10 times greater (Verbunt et al. 2017;Shamohammadi et al. 2024).Therefore we feel that that the larger DM variations observed in the TPA pulsars can be reasonably attributed to the fact that the line of sight to the pulsar traverses the ISM at a rate up to an order of magnitude faster than for the PTA sample.If one assumes that the observed DM slope scales linearly with pulsar velocity, and that the dependence on DM is also close to linear over the typical pulsar distance, we suggest a 'rule of thumb' for predicting the magnitude of DM variations in a pulsar by |DMslope| ∼ 0.01 cm −3 pc yr −1 DM 100 cm −3 pc where  is the pulsar's velocity.We note that in reality there are large, and potentially systematic variations in this value due to the complex structure of the IISM.

RM slopes vs DM
There are 15 pulsars for which we measure a significant DM slope and a significant RM slope.These are plotted in the upper panel of Figure 6.There is a clear correlation between the magnitude of the RM and DM slopes in these pulsars.There does not appear to be any correlation between the sign of the DM slope and the relative sign of the RM slope (i.e. if the magnitude of RM is increasing or decreasing) amongst these pulsars.Put together, this suggests that for these pulsars at least, the DM slopes and RM slopes are likely arising from electron density variations in the same structures in the IISM, though these structures do not dominate the overall RM.The lower panel of Figure 6 also shows that more broadly across the TPA sample, the amplitude of the RM gradients do not scale strongly with DM, perhaps closer to √ DM or even flat, rather than the linear growth seen in the DM gradients.A scaling with √ DM may suggest that the RM slopes originate from a sum over discrete structures in the ISM, perhaps dominated by variations in the local magnetic fields which would be transparent to the DM.A sum over discrete structures could be expected to increase as the square root of the number of structures on the line of sight, and hence roughly as √ DM (Backer et al. 1993).Hence, the observed RM slopes may have two origins, those that arise from the same turbulent IISM effects that cause the DM slopes, and hence scale proportionally with the DM slopes, and those that are caused by discrete magnetised structures on the line of sight and hence scale most closely to the square-root of the distance to the pulsar.

Glitches
Table 2 shows the preliminary results for the glitches in the sample as the fractional change in the spin frequency Δ/ and its error.There are 11 glitches in 9 pulsars.Further studies are ongoing to fully characterise these glitches in terms of transient events or changes in .
Glitches of comparable size to those listed here have previously been seen in 7 of the pulsars (see compilations in Lower et al. 2021 andBasu et al. 2022).For two of the pulsars we have detected glitches for the first time.PSR J1341−6023 is a canonical middleaged pulsar; these pulsars undergo large but rare glitches.In contrast, PSR J1935+2025 has spin parameters similar to the 'Vela-like' pulsars, which undergo quasi-periodic glitches on timescales of years to decades.In our three years of monitoring this pulsar we have detected two glitches, one relatively small in amplitude and the other relatively large.Again, this behaviour is common to the Vela-like pulsars.

Variations in 𝜈
The observed  variations can be broadly classified into 5 classes: (i) timeseries that show a constant slope in , i.e.where there is a significant ; (ii) constant timeseries with no  variations; (iii) timeseries characterized by smooth, long timescale variations; (iv) timeseries that show significant short term variability; and (v) timeseries that show a periodic or quasi-periodic signature.
Four examples of () behaviour are shown in Figure 7. PSR J1347−5947 shows a low-level stochastic wandering of , seemingly typical behaviour for this type of moderate  pulsars.PSR J1531−5610 shows a straight line in , indicative of a large value of .Indeed,  was measured for this pulsar by Parthasarathy et al. (2020) who obtained a value of 1.37(2) × 10 −23 s −3 .The value we measure here, 1.23(5) × 10 −23 s −3 over a much shorter time-span is similar to their value.PSR J1602−5100 was shown by Brook et al. (2016) to undergo a profile change accompanied by a step change in  around MJD 54700.In the MeerKAT data we see a similar  event around MJD 59200 which also lasts for several hundred days.Unlike the event in Brook et al. (2016), here the  returns to its original preevent value.PSR J1638−5226 is an example of a pulsar with highly regular quasi-periodic oscillations in  with a time-scale of some 220 days.Such quasi-periodic behaviour is seen widely across the pulsar population (Hobbs et al. 2010;Niţu et al. 2022) and thought to be linked to emission changes (Lyne et al. 2010).

CONCLUSIONS
This paper presents the first data release from the MeerTime TPA, and we highlight a few of the most striking features from the 597  pulsars.A range of spin-down behaviours is observed, and future work will attempt to identify correlations with profile shape changes (or other properties).We also observe significant DM and RM variations in 15% and 13% of the pulsars over the 4-year timescale of our observations.These DM variations seem larger than a similar sample of MSPs, and we attribute this mainly to the different velocity dispersion in these two families of pulsars, as the line of sight traverses the ISM more quickly in the normal pulsar population.A full analysis of all the data in this data release is outside the scope of this data release paper, but we anticipate that the data will be used for a wide range of projects within the TPA and externally.

Figure 1 .
Figure 1.Period vs period derivative (- diagram) for the 597 pulsars in this data release (blue), the full TPA census (pink) and all pulsars in the ATNF pulsar catalogue (psrcat, version 1.70, https://www.atnf.csiro.au/research/pulsar/psrcat).Note that not all pulsars in psrcat are shown due to the plot axis limits.

)Figure 2 .
Figure 2. Top panel shows the timing residuals as a function of epoch for PSR PSR J1117−6154, followed by the values of , the flux density, the DM, and ionosphere corrected RM.The solid lines show the best fitting slopes to the DM and RM versus time.1- errors are shown on all panels, but too small to be visible on residual and flux density.The large error bar on the DM fit represents an estimate of the systematic error in the mean DM due to frequency evolution of the pulse profile (see Section 3.3 for a description of how this is estimated).

Figure 4 .
Figure 4. DM slope as a function of DM for the MeerTime TPA pulsars, and the MSPs in the MeerTime PTA data release.Downward triangles mark 2- upper limits.Diagonal dashed lines of proportionality are shown for reference.The dot dashed curved line shows the proportionality relationship derived from scattering measurements from Equation 9.

Figure 5 .
Figure 5. DM slope, scaled by DM, as a function of velocity.For TPA pulsars velocity is inferred from proper motion measurements and distance estimates taken from psrcat.PTA velocities taken from Shamohammadi et al. (2024).Not all pulsars have proper motions measured.Downward triangles mark 2- upper limits.Pulsars with parallax derived distances are shown with filled symbols, those with DM-derived distances are shown with open symbols.

Figure 6 .
Figure 6.Upper panel: Magnitude of RM slope vs DM slope for the 15 pulsarswe measure both RM and DM slopes.Lower panel: Magnitude of RM slope vs DM for all pulsars where RM slopes were detected (red points) and those with both RM and DM slopes (blue crosses).In both panels the black gray dashed line is a line of linear proportionality.

Table 1 .
Description of the columns in the observation tables.Header lines give the file extension for each table.Epochs are topocentric values close to the midpoint of the observation.