ISM

    interstellar medium
IISM

    ionised interstellar medium
DM

    dispersion measure
ToA

    time of arrival
MSP

    millisecond pulsar
PTA

    pulsar timing array
EPTA

    European Pulsar Timing Array
PPTA

    Parkes Pulsar Timing Array
NANOGrav

    North American Nanohertz Observatory for Gravitational Waves
IPTA

    International Pulsar Timing Array
LOFAR

    LOw Frequency ARray
ILT

    International LOw Frequency ARray (LOFAR) Telescope
HBA

    high-band antenna
LBA

    low-band antenna
GLOW

    German LOng Wavelength
MPIfR

    Max-Planck Institut für Radioastronomie
LuMP

     LOFAR und Max-Planck Institut für Radioastronomie (MPIfR) Pulsare
RFI

    radio frequency interference
S/N

    signal-to-noise ratio
PulP

    the pulsar pipeline
11institutetext:  Max-Planck-Institut für Radioastronomie, Auf dem Hügel 69,
53121 Bonn, Germany 22institutetext: Fakultät für Physik, Universität
Bielefeld, Postfach 100131, 33501 Bielefeld, Germany 33institutetext: ASTRON,
the Netherlands Institute for Radio Astronomy, Oude Hoogeveensedijk 4, 7991 PD
Dwingeloo, The Netherlands 44institutetext: Gravitational Wave Data Centre,
Swinburne University of Technology, P.O. Box 218, Hawthorn, VIC 3122,
Australia 55institutetext:  Centre for Astrophysics and Supercomputing,
Swinburne University of Technology, P.O. Box 218, Hawthorn, VIC 3122,
Australia 66institutetext:  LPC2E – Université d’Orléans / CNRS, 45071 Orléans
cedex 2, France 77institutetext: Station de Radioastronomie de Nançay,
Observatoire de Paris, PSL Research University, CNRS, Université d’Orléans,
OSUC, 18330 Nançay, France 88institutetext: SKA South Africa, The Park, Park
Road, Pinelands 7405, South Africa 99institutetext: Department of Physics and
Astronomy, University of the Western Cape, Private Bag X17, Bellville 7535,
South Africa 1010institutetext: Jodrell Bank Centre for Astrophysics,
Department of Physics and Astronomy, The University of Manchester, Manchester
M13 9PL, UK 1111institutetext: Technische Universität Berlin, Institut für
Geodäsie und Geoinformationstechnik, Fakultät VI, Sekr. H 12, Straße des 17.
Juni 135, 10623 Berlin, Germany 1212institutetext: GFZ German Research Centre
for Geosciences, Telegrafenberg, 14473 Potsdam, Germany 1313institutetext: SKA
Organisation, Jodrell Bank, Macclesfield SK11 9FT, UK 1414institutetext:
Astro Space Centre, Lebedev Physical Institute, Russian Academy of Sciences,
Profsoyuznaya Str. 84/32, Moscow 117997, Russia 1515institutetext: CSIRO
Astronomy and Space Science, PO Box 1130 Bentley, WA 6102, Australia
1616institutetext: Canadian Institute for Theoretical Astrophysics, University
of Toronto, 60 St. George Street, Toronto, ON M5S 3H8, Canada
1717institutetext: Anton Pannekoek Institute, University of Amsterdam, Postbus
94249, 1090 GE Amsterdam, The Netherlands 1818institutetext:  Hamburger
Sternwarte, University of Hamburg, Gojenbergsweg 112, 21029 Hamburg, Germany
1919institutetext:  Max-Planck-Institut für Astrophysik, Karl-Schwarzschild-
Straße 1, 85748 Garching b. München, Germany 2020institutetext:  Thüringer
Landessternwarte, Sternwarte 5, 07778 Tautenburg, Germany 2121institutetext:
Leibniz-Institut für Astrophysik Potsdam (AIP), An der Sternwarte 16, 14482
Potsdam, Germany

# Dispersion measure variability for 36 millisecond pulsars at 150 MHz with
LOFAR

J. Y. Donner Dispersion measure variability for 36 millisecond pulsars at 150
MHz with LOFARDispersion measure variability for 36 millisecond pulsars at 150
MHz with LOFARDispersion measure variability for 36 millisecond pulsars at 150
MHz with LOFARDispersion measure variability for 36 millisecond pulsars at 150
MHz with LOFAR J. P. W. Verbiest Dispersion measure variability for 36
millisecond pulsars at 150 MHz with LOFARDispersion measure variability for 36
millisecond pulsars at 150 MHz with LOFARDispersion measure variability for 36
millisecond pulsars at 150 MHz with LOFARDispersion measure variability for 36
millisecond pulsars at 150 MHz with LOFAR C. Tiburzi Dispersion measure
variability for 36 millisecond pulsars at 150 MHz with LOFARDispersion measure
variability for 36 millisecond pulsars at 150 MHz with LOFAR S. Osłowski
Dispersion measure variability for 36 millisecond pulsars at 150 MHz with
LOFARDispersion measure variability for 36 millisecond pulsars at 150 MHz with
LOFARDispersion measure variability for 36 millisecond pulsars at 150 MHz with
LOFARDispersion measure variability for 36 millisecond pulsars at 150 MHz with
LOFAR J. Künsemöller Dispersion measure variability for 36 millisecond pulsars
at 150 MHz with LOFARDispersion measure variability for 36 millisecond pulsars
at 150 MHz with LOFAR A.-S. Bak Nielsen Dispersion measure variability for 36
millisecond pulsars at 150 MHz with LOFARDispersion measure variability for 36
millisecond pulsars at 150 MHz with LOFARDispersion measure variability for 36
millisecond pulsars at 150 MHz with LOFARDispersion measure variability for 36
millisecond pulsars at 150 MHz with LOFAR J.-M. Grießmeier Dispersion measure
variability for 36 millisecond pulsars at 150 MHz with LOFARDispersion measure
variability for 36 millisecond pulsars at 150 MHz with LOFARDispersion measure
variability for 36 millisecond pulsars at 150 MHz with LOFARDispersion measure
variability for 36 millisecond pulsars at 150 MHz with LOFAR M. Serylak
Dispersion measure variability for 36 millisecond pulsars at 150 MHz with
LOFARDispersion measure variability for 36 millisecond pulsars at 150 MHz with
LOFARDispersion measure variability for 36 millisecond pulsars at 150 MHz with
LOFARDispersion measure variability for 36 millisecond pulsars at 150 MHz with
LOFAR M. Kramer Dispersion measure variability for 36 millisecond pulsars at
150 MHz with LOFARDispersion measure variability for 36 millisecond pulsars at
150 MHz with LOFARDispersion measure variability for 36 millisecond pulsars at
150 MHz with LOFARDispersion measure variability for 36 millisecond pulsars at
150 MHz with LOFAR J. M. Anderson Dispersion measure variability for 36
millisecond pulsars at 150 MHz with LOFARDispersion measure variability for 36
millisecond pulsars at 150 MHz with LOFARDispersion measure variability for 36
millisecond pulsars at 150 MHz with LOFARDispersion measure variability for 36
millisecond pulsars at 150 MHz with LOFAR O. Wucknitz Dispersion measure
variability for 36 millisecond pulsars at 150 MHz with LOFARDispersion measure
variability for 36 millisecond pulsars at 150 MHz with LOFAR E. Keane
Dispersion measure variability for 36 millisecond pulsars at 150 MHz with
LOFARDispersion measure variability for 36 millisecond pulsars at 150 MHz with
LOFAR V. Kondratiev Dispersion measure variability for 36 millisecond pulsars
at 150 MHz with LOFARDispersion measure variability for 36 millisecond pulsars
at 150 MHz with LOFARDispersion measure variability for 36 millisecond pulsars
at 150 MHz with LOFARDispersion measure variability for 36 millisecond pulsars
at 150 MHz with LOFAR C. Sobey Dispersion measure variability for 36
millisecond pulsars at 150 MHz with LOFARDispersion measure variability for 36
millisecond pulsars at 150 MHz with LOFAR J. W. McKee Dispersion measure
variability for 36 millisecond pulsars at 150 MHz with LOFARDispersion measure
variability for 36 millisecond pulsars at 150 MHz with LOFAR A. V. Bilous
Dispersion measure variability for 36 millisecond pulsars at 150 MHz with
LOFARDispersion measure variability for 36 millisecond pulsars at 150 MHz with
LOFAR R. P. Breton Dispersion measure variability for 36 millisecond pulsars
at 150 MHz with LOFARDispersion measure variability for 36 millisecond pulsars
at 150 MHz with LOFAR M. Brüggen Dispersion measure variability for 36
millisecond pulsars at 150 MHz with LOFARDispersion measure variability for 36
millisecond pulsars at 150 MHz with LOFAR B. Ciardi Dispersion measure
variability for 36 millisecond pulsars at 150 MHz with LOFARDispersion measure
variability for 36 millisecond pulsars at 150 MHz with LOFAR M. Hoeft
Dispersion measure variability for 36 millisecond pulsars at 150 MHz with
LOFARDispersion measure variability for 36 millisecond pulsars at 150 MHz with
LOFAR J. van Leeuwen Dispersion measure variability for 36 millisecond pulsars
at 150 MHz with LOFARDispersion measure variability for 36 millisecond pulsars
at 150 MHz with LOFARDispersion measure variability for 36 millisecond pulsars
at 150 MHz with LOFARDispersion measure variability for 36 millisecond pulsars
at 150 MHz with LOFAR C. Vocks Dispersion measure variability for 36
millisecond pulsars at 150 MHz with LOFARDispersion measure variability for 36
millisecond pulsars at 150 MHz with LOFAR

(Received 24 Sep 2020 / Accepted 15 Nov 2020)

###### Abstract

Context.  Radio pulses from pulsars are affected by plasma dispersion, which
results in a frequency-dependent propagation delay. Variations in the
magnitude of this effect lead to an additional source of red noise in pulsar
timing experiments, including pulsar timing arrays that aim to detect
nanohertz gravitational waves.

Aims.  We aim to quantify the time-variable dispersion with much improved
precision and characterise the spectrum of these variations.

Methods.  We use the pulsar timing technique to obtain highly precise
dispersion measure (DM) time series. Our dataset consists of observations of
36 millisecond pulsars, which were observed for up to 7.1 years with the LOFAR
telescope at a centre frequency of $\sim 150$ MHz. Seventeen of these sources
were observed with a weekly cadence, while the rest were observed at monthly
cadence.

Results.  We achieve a median DM precision of the order of $10^{-5}$ cm-3 pc
for a significant fraction of our sources. We detect significant variations of
the DM in all pulsars with a median DM uncertainty of less than $2\times
10^{-4}$ cm-3 pc. The noise contribution to pulsar timing experiments at
higher frequencies is calculated to be at a level of
0.1–$10\text{\,}\mathrm{\SIUnitSymbolMicro s}$ at 1.4 GHz over a timespan of a
few years, which is in many cases larger than the typical timing precision of
$1\text{\,}\mathrm{\SIUnitSymbolMicro s}$ or better that PTAs aim for. We
found no evidence for a dependence of DM on radio frequency for any of the
sources in our sample.

Conclusions.  The DM time series we obtained using LOFAR could in principle be
used to correct higher-frequency data for the variations of the dispersive
delay. However, there is currently the practical restriction that pulsars tend
to provide either highly precise times of arrival at 1.4 GHz or a high DM
precision at low frequencies, but not both, due to spectral properties.
Combining the higher-frequency ToAs with those from LOFAR to measure the
infinite-frequency ToA and DM would improve the result.

###### Key Words.:

pulsars: general – ISM: structure – gravitational waves

##  1 Introduction

Pulsars are highly magnetised, rapidly rotating neutron stars, the remnants of
massive stars that ended their lives in a supernova explosion. Pulsars emit
beams of electromagnetic radiation, pronounced mostly at radio frequencies,
which sweep around in space as the neutron star rotates. If any of the beams
cross the line of sight to the Earth, we detect regular pulses of radiation.
While the shape of individual pulses can vary significantly, the average pulse
shape, integrated over hundreds of pulses, is usually stable, which allows for
a precise measurement of the ToAs of the integrated pulses.

There are two major distinct populations of pulsars: the canonical pulsars
with pulse periods of the order of $\sim$0.1 s to $\sim$20 s and the MSPs,
with typical spin periods of 1.4 ms to $\sim$30 ms. The latter were ‘spun up’
to shorter pulse periods by the transfer of mass and angular momentum from a
binary stellar companion, which is why they are also called ‘recycled’ pulsars
(Tauris et al. 2015). Millisecond pulsars are of particular scientific
interest as the much shorter pulse periods allow for a more precise
determination of the ToAs, and their rotation was found to be much more stable
than that of canonical pulsars (Hobbs et al. 2010; Verbiest et al. 2009).

#### Pulsar timing.

The high rotational stability of pulsars, and MSPs in particular, together
with the precise measurements afforded by pulsar timing, allows for extremely
accurate modelling of the pulsar’s astrometric and astrophysical properties
with a so-called timing model, a method called pulsar timing (see Lorimer &
Kramer 2005). The model is usually evaluated by inspecting the timing
residuals, that is, the difference between the measured ToAs of the pulses and
the ones predicted by the model.

#### Pulsar timing arrays.

A major application of pulsars is in the so-called PTA projects, which aim to
detect nanohertz frequency gravitational waves (see, e.g. Hobbs & Dai 2017;
Tiburzi 2018; Burke-Spolaor et al. 2019). The basic idea behind these
experiments is that passing gravitational waves would distort the spacetime
around Earth in a way that the arrival times of some pulsars would be delayed,
while the pulses of other pulsars arrive earlier at the telescope. By
observing a large number of spatially separated, precisely timed pulsars over
a time span of years to decades, a gravitational wave would be measurable in
the spatial correlation of the timing residuals of different pulsars. One
major focus of the PTAs is to detect a stochastic gravitational wave
background signature such as the Hellings and Downs curve (Hellings & Downs
1983). This background is expected to be caused by inspiralling super-massive
black hole binaries. To achieve a high precision in these experiments, many
sources of noise have to be taken into account, including propagation effects
in the ionised interstellar medium (IISM).

#### Influence of the interstellar medium.

The IISM is a cold plasma of ionised particles. An electromagnetic wave
passing through this medium will experience a frequency-dependent change in
group velocity, a phenomenon called _dispersion_. Specifically, this induces
an additional delay $\Delta t$ in the ToAs, which can be approximated as (see
Lorimer & Kramer 2005):

$\Delta t=\frac{\rm DM}{k\nu^{2}},$ (1)

where $k\simeq 2.41\times 10^{-4}$ cm-3 pc MHz-2 s-1 is the dispersion
constant111Although $k$ can be calculated more precisely, we use this rounded
value as it is common practice in the field. See Kulkarni (2020) for a
discussion., $\nu$ the observing frequency (expressed in MHz) and DM is the
‘dispersion measure’ (the integrated electron column density along the line of
sight, expressed in cm-3 pc), which is defined as:

${\rm DM}=\int_{0}^{d}n_{\rm e}(l){\rm d}l,$ (2)

where $d$ is the distance to the pulsar and $n_{\rm e}$ is the electron
density.

The IISM is inhomogeneous, and so as our line of sight to the pulsar changes
due to its proper motion, the Earth’s motion and the motion of the IISM we
sample different regions of the IISM, which implies a temporal variability in
the DM. The variations in the DM lead to a variable dispersive delay (see Eq.
1) and, thus, are a major source of noise in pulsar timing.

Using relatively high observing frequencies (¿ 1.4 GHz) minimises the impact
of DM variations on the ToAs as the dispersive delay scales with $\nu^{-2}$
(see Eq. 1). However, the IISM turbulence spectrum is steep with significantly
more power at larger scales (Armstrong et al. 1995), so for long timing
campaigns, the dispersive delays will sooner or later still have a significant
impact on the ToAs. Also, the observing frequency cannot be increased
indefinitely, as pulsars have rather steep spectra (Maron et al. 2000), which
leads to a loss of signal-to-noise ratio (S/N) at very high frequencies. Out
of the 33 MSPs in their sample, Lazarus et al. (2016) selected 12 candidates
with promising spectral indices to run test observations at 5 GHz (3 at 9
GHz). All the selected sources are detected, but the S/N is significantly
lower than at 1.4 GHz, especially in the 9 GHz band. However, some of the
sources have a sub-$\mathrm{\SIUnitSymbolMicro s}$ timing precision (post-fit
RMS) at 5 GHz, so this approach can work for some flat-spectrum sources.

A more generally applicable solution is to correct for the time-variable
dispersive delays by subtracting them from the ToAs, which effectively yields
ToAs at infinite frequency. Dispersive delays can be measured either by using
multi-frequency observations with several receivers, or by splitting wide-band
observations into multiple sub-bands, but this kind of observation is not
always available in PTA experiments (see, e.g. Desvignes et al. 2016). Also,
the correction for the DM increases the uncertainty of the infinite-frequency
ToA and this effect is especially large at high frequencies and small
fractional bandwidths.

Another approach is to use low-frequency observations to measure the DM very
precisely and use these measurements to correct the higher-frequency ToAs,
which are often more sensitive and less affected by IISM effects. One
potential complication of this approach is the possibility of frequency-
dependent DMs (Cordes et al. 2016; Donner et al. 2019), which are caused by
the fact that due to interstellar scattering, low-frequency observations
effectively sample a larger volume of the IISM and are sensitive to IISM
structures on a larger scale than at high frequencies. Therefore, dispersive
delays measured at low frequency may not be representative for those
experienced at high frequencies. Additionally, time-variable scattering can
cause variations in the shape of the pulse profile that, if not corrected for,
cause a frequency-dependent delay in the timing. While scattering induces a
delay scaling with $\nu^{-4}$, some of its signature can be absorbed into a DM
measurement, so this effect can cause a mismatch of the apparent DM at
different frequencies.

In this paper, we present precise DM time series of 36 MSPs using low-
frequency ($\sim$150 MHz) observations with a particular focus on pulsars used
in PTA experiments. In Sect. 2 we describe our observational setup and our
sources, while Sect. 3 explains our data analysis. We discuss our findings in
Sect. 4 and present our conclusions in Sect. 5.

##  2 Observations

We used observations taken with the International  Telescope (ILT) high-band
antennas, described in detail by van Haarlem et al. (2013). Specifically, we
used the data from two different pulsar monitoring campaigns, both of which
are still ongoing. Campaign 1 uses the LOFAR Core stations situated in the
north-east of the Netherlands to observe the sources of interest with monthly
cadence for up to 7.1 years between 19 December 2012 and 14 January 2020.
These data were coherently dedispersed and reduced using the pulsar pipeline
(PulP) described by Stappers et al. (2011) and Kondratiev et al. (2016). It
produces data cubes with resolution in frequency (195.3125 kHz-wide channels),
time (10-sec sub-integrations), polarisation (four coherency products) and
rotational phase (256 to 1024 phase bins). The data are stored in timer
format, which is similar to the psrfits format described in Hotan et al.
(2004). Each observation uses a centre frequency of 148.9 MHz, a bandwidth of
78.1 MHz (400 frequency channels) and lasts 5 to 30 minutes, depending on the
brightness of the pulsars.

In Campaign 2, the six LOFAR stations of the German LOng Wavelength (GLOW)
consortium, located in Effelsberg (telescope identifier DE601),
Unterweilenbach (DE602), Tautenburg (DE603), Potsdam-Bornim (DE604), Jülich
(DE605) and Norderstedt (DE609), were used as individual stand-alone
telescopes, not connected to the ILT network. The beamformed data from the
stations were sent to the MPIfR and the Forschungszentrum Jülich on dedicated
high-speed links, where recording computers ran the dedicated  und  Pulsare
(LuMP)222Publicly available at https://github.com/AHorneffer/lump-lofar-und-
mpifr-pulsare and described on https://deki.mpifr-
bonn.mpg.de/Cooperations/LOFAR/Software/LuMP. data-taking software. LuMP
formats and otherwise prepares the beamformed pulsar data for subsequent (off-
line but near-real-time) phase-resolved averaging (commonly referred to as
‘folding’) using the dspsr software package (van Straten & Bailes 2011). This
produces data cubes in the same format as PulP. The data we used from this
campaign were taken between 20 August 2013 and 08 January 2020. The resulting
dataset of these observations covers a time span of up to 6.4 years per
pulsar, with a weekly cadence and typical integration times of 1 to 3 hours.
Due to the difference in collecting areas (a factor of $\sim$10), the faintest
pulsars were only observed with the Core and are not detectable with the
international stations. Early observations and the observations of DE601 have
a total bandwidth of 95.3 MHz (488 frequency channels), centred at 149.9 MHz.
For technical reasons, the bandwidth of the other stations was reduced to 71.5
MHz (366 channels) in February 2015. In order to minimise the impact of the
bandwidth reduction on the scientific quality of our data, the centre
frequency was shifted to align the observed bandwidth with the most sensitive
part of the bandpass, resulting in a new centre frequency of 153.8 MHz. This
implies a shift in centre frequency by an integer number of frequency channels
(20), so that the frequencies of individual channels remained constant over
the entire dataset.

Table 1 shows detailed information on the observation characteristics for each
source. We excluded PSR J1939+2134 (PSR B1937+21) from our analysis because
its strongly variable scattering has a significant impact on the DM estimation
(e.g. it shows the largest variation in scattering time of all sources in
Levin et al. 2016). Also, due to the very strong scattering, the scattering
tail of the pulse profile merges with the next pulse, so the profile is very
wide (see Fig. 2 of Kondratiev et al. 2016). Observations of this pulsar at
higher frequencies or more advanced analysis techniques like cyclic
spectroscopy (Demorest 2011) could lead to more robust and similarly precise
DM measurements.

Table 1: Summary of observations. Given are the pulsar name in J2000
coordinates, the pulse period $P$, the catalogue DM, the ecliptic latitude
$\beta$, the total time span of the observations, the number of observations
with GLOW and the LOFAR Core $N_{\text{glow}}$ and $N_{\text{core}}$, the
median frequency-integrated ToA uncertainty, the median uncertainty of
individual DM measurements, the median reduced $\chi^{2}$ of the individual DM
fits, and the number of frequency channels in each observation. The last
column shows whether the pulsar is used in different PTA projects (i.e. E:
European Pulsar Timing Array, EPTA; P: Parkes Pulsar Timing Array, PPTA; N:
North American Nanohertz Observatory for Gravitational Waves, NANOGrav; see
Verbiest et al. 2016; Arzoumanian et al. 2018). A pulsar observed by any of
the PTAs is also used by the International Pulsar Timing Array (IPTA).

source name | $P$ | DM | $\beta$ | span | $N_{\text{glow}}$ | $N_{\text{core}}$ | med($\sigma_{\text{ToA}}$) | med($\sigma_{\text{DM}}$) | med$\left(\frac{\chi^{2}}{n_{\text{free}}}\right)$ | $N_{\text{chan}}$ | PTA
---|---|---|---|---|---|---|---|---|---|---|---
(J2000) | (ms) | $\left(\frac{\text{pc}}{\text{cm}^{3}}\right)$ | (deg) | (yrs) |  |  | ($\mathrm{\SIUnitSymbolMicro s}$) | $\left(10^{-5}\frac{\text{pc}}{\text{cm}^{3}}\right)$ |  |  |
J0030+0451 | $4.9$ | $4.3$ | $1.4$ | $7.0$ | 355 | 79 | 7 | 16 | $1.0$ | 10 | E P N
J0034$-$0534 | $1.9$ | $13.8$ | $-8.5$ | $7.0$ | 1041 | 78 | 1 | 3 | $1.0$ | 20 | E P N
J0218+4232 | $2.3$ | $61.2$ | $27.0$ | $6.5$ | 539 | 74 | 13 | 35 | $0.9$ | 10 | E P N
J0407+1607 | $25.7$ | $35.6$ | $-4.7$ | $6.4$ | 386 | 65 | 24 | 56 | $1.0$ | 10 | -
J0621+1002 | $28.9$ | $36.5$ | $-13.3$ | $7.0$ | 10 | 74 | 251 | 568 | $1.3$ | 10* | E P N
J0645+5158 | $8.9$ | $18.2$ | $28.9$ | $6.9$ | 268 | 4 | 3 | 7 | $0.9$ | 10 |  E P N
J0740+6620 | $2.9$ | $15.0$ | $44.1$ | $4.8$ | 261 | 0 | 2 | 4 | $0.9$ | 10 |  E P N
J0751+1807 | $3.5$ | $30.2$ | $-2.8$ | $7.1$ | 0 | 71 | 34 | 89 | $0.6$ | 10* | E P N
J0952$-$0607 | $1.4$ | $22.4$ | $-17.9$ | $2.9$ | 37 | 33 | 9 | 14 | $1.0$ | 8* | -
J1012+5307 | $5.3$ | $9.0$ | $38.8$ | $7.1$ | 1210 | 79 | 4 | 12 | $1.2$ | 10 | E P N
J1022+1001 | $16.5$ | $10.3$ | $-0.1$ | $7.1$ | 1224 | 80 | 5 | 12 | $1.2$ | 20 | E P N
J1024$-$0719 | $5.2$ | $6.5$ | $-16.0$ | $7.1$ | 5 | 81 | 32 | 80 | $0.9$ | 8* | E P N
J1125+7819 | $4.2$ | $11.2$ | $62.5$ | $4.8$ | 282 | 0 | 3 | 8 | $0.9$ | 10 |  E P N
aJ1300+1240a | $6.2$ | $10.2$ | $17.6$ | $7.1$ | 377 | 80 | 3 | 7 | $0.8$ | 20 | -
J1400$-$1431 | $3.1$ | $4.9$ | $-2.2$ | $4.3$ | 145 | 49 | 7 | 6 | $0.8$ | 8* | -
J1544+4937 | $2.2$ | $23.2$ | $65.9$ | $5.6$ | 5 | 53 | 19 | 59 | $1.0$ | 6* | -
J1552+5437 | $2.4$ | $22.9$ | $70.7$ | $3.7$ | 0 | 41 | 3 | 7 | $0.9$ | 10* | -
J1640+2224 | $3.2$ | $18.4$ | $44.1$ | $7.1$ | 435 | 79 | 7 | 18 | $0.8$ | 8 | E P N
J1658+3630 | $33.0$ | $3.0$ | $58.7$ | $2.6$ | 164 | 6 | 7 | 18 | $0.9$ | 10 | -
J1713+0747 | $4.6$ | $16.0$ | $30.7$ | $7.0$ | 7 | 79 | 12 | 29 | $0.6$ | 10* | E P N
J1730$-$2304 | $8.1$ | $9.6$ | $0.2$ | $6.5$ | 0 | 74 | 21 | 49 | $0.8$ | 10* | E P N
J1738+0333 | $5.9$ | $33.8$ | $26.9$ | $7.1$ | 0 | 77 | 39 | 110 | $0.7$ | 6* | E P N
J1744$-$1134 | $4.1$ | $3.1$ | $11.8$ | $7.0$ | 372 | 79 | 14 | 19 | $0.9$ | 10 | E P N
J1853+1303 | $4.1$ | $30.6$ | $35.7$ | $6.5$ | 1 | 72 | 21 | 50 | $0.8$ | 8* | E P N
bJ1857+0943b | $5.4$ | $13.3$ | $32.3$ | $7.1$ | 0 | 75 | 41 | 82 | $1.3$ | 10* | E P N
J1911$-$1114 | $3.6$ | $31.0$ | $11.1$ | $7.0$ | 0 | 74 | 32 | 68 | $0.9$ | 10 | E P N
J1918$-$0642 | $7.6$ | $26.6$ | $15.4$ | $7.0$ | 0 | 75 | 55 | 136 | $0.8$ | 10* | E P N
J1923+2515 | $3.8$ | $18.9$ | $46.7$ | $6.6$ | 0 | 73 | 6 | 15 | $0.9$ | 10* |  E P N
J1944+0907 | $5.2$ | $24.4$ | $29.9$ | $6.5$ | 35 | 64 | 33 | 80 | $0.8$ | 6* |  E P N
cJ1955+2908c | $6.1$ | $104.5$ | $48.7$ | $6.6$ | 0 | 74 | 100 | 273 | $0.9$ | 10* | E P N
J2043+1711 | $2.4$ | $20.7$ | $34.0$ | $6.7$ | 3 | 71 | 6 | 17 | $0.9$ | 10* |  E P N
J2051$-$0827 | $4.5$ | $20.7$ | $8.8$ | $6.5$ | 6 | 74 | 14 | 29 | $0.7$ | 10* | -
J2145$-$0750 | $16.1$ | $9.0$ | $5.3$ | $7.0$ | 1010 | 81 | 3 | 9 | $1.1$ | 20 | E P N
J2222$-$0137 | $32.8$ | $3.3$ | $8.0$ | $3.8$ | 130 | 41 | 22 | 49 | $0.9$ | 10* | -
J2302+4442 | $5.2$ | $13.7$ | $45.7$ | $6.2$ | 0 | 71 | 20 | 51 | $0.8$ | 10* |  E P N
J2317+1439 | $3.4$ | $21.9$ | $17.7$ | $7.0$ | 381 | 77 | 2 | 6 | $1.0$ | 10 | E P N

* For these pulsars, an analytic standard template without frequency resolution was used.
a PSR B1257+12
b PSR B1855+09
c PSR B1953+29

##  3 Data analysis

###  3.1 Pre-processing

The basic processing in this work has been carried out with the psrchive
(Hotan et al. 2004; van Straten et al. 2012) software package. As a first
step, the data were cleaned from radio frequency interference (RFI), by using
a modified version of the ‘surgical’ algorithm of the clean.py script from the
coastguard (Lazarus et al. 2016) python package333The version we used was
provided by Künkel (2017) and is publicly available at
https://github.com/larskuenkel/iterative_cleaner. to remove affected frequency
channels and sub-integrations. In the rare case of outliers due to remaining
RFI, the observations were also manually inspected and additional cleaning was
applied using the psrchive program pazi. On average, 9.0% of data were removed
in this process.

Pulse profiles can be polarisation dependent, so the total intensity profile
can be significantly distorted if the different polarisation channels are not
calibrated correctly, which affects the ToAs. As LOFAR antennas are static,
there are several time-dependent projection effects to take into account. The
data were calibrated in polarisation following the methods outlined in Noutsos
et al. (2015) using the dreambeam444Publicly available at
https://github.com/2baOrNot2ba/dreamBeam, by T. Carozzi. python package to
calculate the Jones matrices. For some pulsars (especially PSRs J1022+1001 and
J2145$-$0750) this significantly improved the reduced $\chi^{2}$ values of the
DM fits.

To make the dataset as homogeneous as possible, the bandwidth of all
observations has been reduced to 70.3 MHz (i.e. 360 channels, centred at 153.2
MHz). To achieve this, two empty dummy channels had to be added at the top of
the band for the LOFAR Core observations.

###  3.2 Timing

For the 20 pulsars also analysed by the International Pulsar Timing Array
(IPTA), we took the timing model from the first IPTA data release (Verbiest et
al. 2016, data combination B)555Publicly available at http://www.ipta4gw.org/
and removed any existing DM models (including Solar wind) or FD parameters
(which describe a frequency-dependent, non-dispersive delay in the timing
residuals, see Arzoumanian et al. 2015), as we are interested in the time
evolution of the DM and using frequency-resolved timing removes the need for
FD parameters. To account for the different time span over which the IPTA
ephemerides were derived, and because the DM estimates are expected to be more
precise at the low radio frequencies of LOFAR, we performed an initial timing
analysis over our dataset, with the aim of updating the reference DM for each
pulsar. For this initial timing analysis, we derived an analytic template from
a single, bright observation by fitting a series of von Mises functions (see,
e.g. Jammalamadaka & SenGupta 2001) to the total intensity profile of the
observation using the program paas. The functions are implemented in the form:

$f(x)=A\cdot e^{\kappa\cdot(\cos(2\pi(x-\mu))-1)},$ (3)

with $A$ being the amplitude of the component, $\kappa$ the compactness and
$\mu$ the pulse phase. The phase $x$ is defined such that one full rotation
corresponds to $x=1$. To obtain the ToAs, we used the fdm666This algorithm is
identical to that described by Taylor (1992), except for the uncertainties.
fdm uses either formal uncertainties or Monte-Carlo simulations. We used the
uncertainties determined from Monte-Carlo simulations. algorithm as
implemented in the program pat on copies of our observations with reduced
frequency resolution (integrated down to 10 frequency channels). Using these
initial ToAs, we calculated the DM and its running average for each
observation as described in Sect. 3.3, to be able to phase-align observations
affected by DM variations.

For the 16 pulsars without a published IPTA timing model, we took our initial
timing model from the ATNF Pulsar Catalogue
(psrcat)777http://www.atnf.csiro.au/people/pulsar/psrcat/ by Manchester et al.
(2005). For PSR J1658+3630, there was no timing model available in the
catalogue, so we used the timing model from Sanidas et al. (2019). As these
timing models often did not phase-align our dataset, we used the tempo2
software (Hobbs et al. 2006) to fit for parameters describing the pulsar’s
rotation and, if needed, the position of the pulsar and its orbit in binary
systems.

Using the optimised timing model and initial DM time series, we phase-aligned
our dataset and created a standard pulse profile template for each pulsar by
averaging only the observations from the observing site with the most
observations of that pulsar. As the S/N of most observations is rather low
(rarely above 100), we integrated the template in frequency to between 6 and
20 channels, depending on the pulsar’s S/N. We then applied a wavelet
smoothing algorithm to the template with the program psrsmooth (Demorest et
al. 2013), to avoid self-standarding (see Appendix A of Hotan et al. 2005).
With this template, we again used the fdm algorithm to calculate the ToAs,
matching each frequency channel of the observation (that was frequency-
integrated accordingly) with the corresponding channel of the template. With
this method, any constant frequency dependence of the pulse profile is
modelled by the template and will not affect our measurement of DM.

If the templates generated in the procedure above were too noisy (due to low
intrinsic flux of the pulsar and a small number of observations), the
smoothing produced unphysical features in the profile shape. In these cases we
used an analytic template without frequency resolution, created in the same
way as in the initial analysis. In principle, this would imply the need of FD
parameters in the timing model to account for the frequency-dependent profile
shape, but in practice the frequency dependence is not significant,
demonstrated by the median $\chi^{2}$ of the DM fits being close to unity (see
Table 1).

###  3.3 Calculation of the DM

To get a precise time series of the DM, we used tempo2 to fit for DM for each
observation individually. This approach avoids any correlation of the measured
DMs with other (time-dependent) timing parameters. To mitigate the impact of
outlier ToAs on our measured DM, we apply an automatic ToA rejection scheme as
was done similarly in Tiburzi et al. (2019). We reject a ToA if its residual-
to-uncertainty ratio $k$ is larger than four times the mean $k$. Specifically,
we repeatedly applied the following rules until no ToA was removed in the
iteration:

We calculated the timing residuals $r_{i}$ for each channel $i$ of a given
observation using the tempo2 general2 plugin and noted the corresponding ToA
uncertainties $\sigma_{i}$. The mean residual-to-uncertainty ratio is then
calculated as $k=\text{mean}(|r_{i}|/\sigma_{i})$ and set to unity if the
result was less than unity. If the condition $|r_{i}|/\sigma_{i}>4k$ was met
for any of the residuals, we rejected the corresponding ToA. Finally, we fit
for DM with tempo2.

The iterative process is necessary as the variable DM gives rise to
significant structure in our frequency-resolved residuals, meaning that many
are far from zero. This is also why the factor $k$ was introduced, as without
it, ToAs at the edges of the band would be removed in the presence of a
dispersive slope. The outliers were caused by a low S/N in the corresponding
frequency channel, which can occur due to the pulsar being intrinsically faint
at that frequency, short integration times, removal of large parts of that
frequency range due to RFI, remaining low-level RFI or a combination of these.
On average, 0.7% of the ToAs of an observation were removed during this step,
mostly at the edges of the band where the telescope is less sensitive. We
estimate that our outlier rejection procedure will produce fewer than ten
false positives across our entire 36-pulsar dataset. While the median reduced
$\chi^{2}$ of the final DM fits is close to unity for all pulsars in our
sample (see Table 1) and the distribution is strongly peaked around 1, there
are observations not following this distribution with reduced $\chi^{2}\gg 1$.
These cases occured when most ToAs of an observation were flawed, usually
because the pulsar was too faint or the entire observation was dominated by
RFI. We excluded all DM measurements with a reduced $\chi^{2}>5$ from our
subsequent analysis (on average 4% of all observations). Finally, we applied
the standard tempo2 procedure of multiplying the DM uncertainty by the square
root of the reduced $\chi^{2}$, because a high reduced $\chi^{2}$ indicates
unmodelled structure or underestimated ToA uncertainties. To be conservative
in our uncertainties and due to the low number of ToAs per fit, we only
applied this procedure for fits with reduced $\chi^{2}>1$.

The resulting DM time series are shown in Figs. 1-3, with the reference DM of
the standard template subtracted. The median DM precision for each source is
shown in Table 1. As the median reduced $\chi^{2}$ is close to unity for all
pulsars in our sample, the data are well fit by the model. From that we can
rule out any significant frequency-dependent structure in the residuals or a
frequency dependence of the DM as discussed in Donner et al. (2019), where we
present a system with more extreme DM variations.

To improve our sensitivity to low-amplitude variations, we computed a running
average of the DM time series. For each observation MJD, a weighted average
over all observations was formed, weighting each DM value by the inverse of
its variance. Additionally, the weights were reduced exponentially with time,
scaled to $1/e$ of their original value over half the averaging window. We
chose an averaging window of 30 days for pulsars with sharp features in their
DM signal (i.e. PSRs J0030+0451, J0034$-$0534, J1022+1001, J1400$-$1431,
J2145$-$0750, and J2317+1439) and 60 days for the rest.

Figure 1: DM time series for 12 MSPs that show a detectable Solar wind signal.
Individual DM measurements are plotted in grey, while running averages and
their uncertainty are represented by solid grey and black lines, respectively.
For faint pulsars only observed with the LOFAR core, the running average is
not plotted due to the low observing cadence. Only points with uncertainties
less than three times the median uncertainty are plotted. The vertical dashed
grey lines indicate the epochs of minimum angular separation from the Sun.
Figure 2: DM time series for six MSPs that were regularly observed with GLOW
and do not show a clear Solar wind signal. The representation is the same as
in Fig. 1.  Figure 3: DM time series for 18 faint MSPs that are only monitored
regularly with the LOFAR Core and do not show a clear Solar wind signal. The
representation is the same as in Fig. 1.

##  4 Discussion

Most of the pulsars in our sample show DM variability, on diverse timescales.
The amplitude of these variations ranges from $\sim 10^{-4}$ cm-3 pc to
$10^{-3}$ cm-3 pc over several years or less. Pulsars that do not show DM
variations (e.g. PSR J1853+1303) are very faint at LOFAR frequencies, with a
single-observation DM uncertainty of $2\times 10^{-4}$ cm-3 pc or worse, so
the non-detection of variations may be caused by a lack of sensitivity, and
variations of the order of a few $10^{-4}$ cm-3 pc cannot be ruled out.

###  4.1 Solar wind

We detect the DM variations due to the Solar wind in 12 pulsars, namely PSRs
J0030+0451, J0034$-$0534, J0407+1607, J0645+5158, J1012+5307, J1022+1001,
J1300+1240, J1400$-$1438, J1730$-$2304, J1744$-$1134, J2145$-$0750, and
J2317+1439 (see Fig. 1). The angular distance of these sources from the
ecliptic ranges from 0.1° (PSR J1022+1001) to 38.8° (PSR J1012+5307, see Table
1). We note that for PSR J1730$-$2304, there are only a few observations close
to the Sun due to the low observing cadence with the LOFAR Core and thus those
observations clearly stand out from the rest of the DM measurements. For some
pulsars with low (¡$10\degree$) absolute ecliptic latitude, we do not clearly
see the impact of the Solar wind due to a lack of measurement precision: PSRs
J0751+1807, J2051$-$0827, and J2222$-$0137.

As shown in Tiburzi et al. (2019), the widely used Solar wind models are
insufficient to correctly account for the highly variable Solar wind, which
leads to short-term variations in the residual DMs of the order of a few
$10^{-4}$ cm-3 pc, even at larger separations from the Sun (up to $\sim
50\degree$). From our study, we can confirm that a spherical model with
constant amplitude (as is usually applied) does not suffice in most cases.
This becomes clear in Fig. 1, where the amplitude of the Solar wind strongly
varies from year to year, which is expected due to the variable Solar
activity. Further analysis of the Solar wind impact on these data will be
presented by Tiburzi et al. (in prep).

###  4.2 Observed DM variations

Lam et al. (2018) report a DM event in PSR J1713+0747 around MJD 57510, where
the DM temporarily drops by about $4\times 10^{-4}$ cm-3 pc. There is no
indication of this event in our LOFAR DM time series (see Fig. 3). This can
partly be explained by a lack of sensitivity as our DM precision is similar to
the amplitude of the event. Additionally, the event is observed at higher
frequencies and could be partly smeared at lower frequencies due to a
frequency-dependent DM, as discussed in Cordes et al. (2016).

Some sources show very linear gradients in DM, for example PSRs J0740+6620,
J1552+5437, and J1400$-$1431, the last of which is also affected by the Solar
wind. While these linear trends are unlikely to frequently occur solely due to
a turbulent process, they can be explained by the pulsar’s motion along the
line of sight (Lam et al. 2016). To quantify overall DM trends, we used a
least-squares fitting routine to fit a linear gradient to the DM time series.
For the sources affected by the Solar wind, we first subtracted a spherically
symmetric Solar wind model with variable amplitude, which was fit year by year
to the difference between the measured DMs and a cubic spline derived from
measurements furthest from the Sun. The resulting DM gradients are given in
Table 2.

Hobbs et al. (2004) found a correlation between DM and its time derivative in
their sample of 374 pulsars. When fitting a power law, they found an exponent
of 0.57(9), consistent with a square-root dependence:
$|\text{dDM/d}t|=0.0002\sqrt{\text{DM}}$ (DM in cm-3 pc and dDM/d$t$ in cm-3
pc yr-1). Applying this analysis to our data, we get an exponent of 0.7(6),
which is consistent with the Hobbs et al. (2004) result, albeit not very
constraining. When we fit the amplitude of a square-root dependence, however,
we get $|\text{dDM/d}t|=0.000024(5)\sqrt{\text{DM}}$, highly inconsistent with
the value from Hobbs et al. (2004), which is an order of magnitude larger.
This suggests that this kind of analysis is biased. One major difference
between the two studies is the sensitivity to DM: Our single-observation
measurements of DM are typically two orders of magnitude more precise than the
overall DM of Hobbs et al. (2004), so we investigated possible correlations
between DM, —dDM/d$t$—, and $\sigma_{\text{DM}}$ in both datasets by fitting
power law relations (see Fig. 4). Notably, the correlation between —dDM/d$t$—
and $\sigma_{\text{DM}}$ in the Hobbs et al. (2004) dataset is nearly as
strong as the correlation between —dDM/d$t$— and DM. Furthermore, only 21 out
of the 374 sources show a —dDM/d$t$— measurement above $3\sigma$ significance,
while two-thirds are below $1\sigma$, so a majority of their measurements are
upper limits that can bias the fit.

Figure 4: Correlations between DM, —dDM/d$t$— and $\sigma_{\text{DM}}$ in our
dataset (panels a to c) and that of Hobbs et al. (2004) (panels d to f). For
each panel, we fit a power law of the form $y=ax^{b}$ to the data and the
slope $b$ of this fit is given in each panel, along with its uncertainty. For
the correlation between DM and $\sigma_{\text{DM}}$, we fit in logspace due to
the lack of an uncertainty in $y$. As it was done by Hobbs et al. (2004), we
excluded all pulsars from their sample with —dDM/d$t$— ¿ 0.01 cm-3 pc yr-1.

Another potential bias lies in the source selection. The sample observed by
Hobbs et al. (2004) is strongly focused on sources in the direction of the
Galactic Centre. As the sensitivity of LOFAR is strongly dependent on source
elevation (see Fig. 1 of Noutsos et al. 2015), the detectability of pulsars
(and therefore our sample) is biased towards high-declination sources.

Table 2: Quantities derived from the DM time series of each pulsar. Given are the average DM derivative over the entire dataset and the value of the structure function $D_{\text{DM}}$ at a time lag of 1000 days. source name | $\left<\frac{\text{dDM}}{\text{d}t}\right>$ | $D_{\text{DM}}(1000\,\text{d})$
---|---|---
(J2000) | $\left(10^{-4}\frac{\text{pc}}{\text{cm}^{3}\,\text{yr}}\right)$ | $\left(10^{-6}\frac{\text{pc}^{2}}{\text{cm}^{6}}\right)$
J0030+0451 | $0.40(7)$ | $0.023$
J0034$-$0534 | $-0.37(2)$ | $0.016$
J0218+4232 | $-6.10(19)$ | $6.2$
J0407+1607 | $-1.30(18)$ | $0.56$
J0621+1002 | 51(4) | $220$
J0645+5158 | $0.02(4)$ | $0.021$
J0740+6620 | $4.40(5)$ | $0.83$
J0751+1807 | $-5.2(6)$ | $2.1$
J0952$-$0607 | $-0.6(2)$ | $0.17$
J1012+5307 | $1.30(3)$ | $0.086$
J1022+1001 | $0.02(5)$ | $0.0038$
J1024$-$0719 | $-1.7(4)$ | $0.90$
J1125+7819 | $1.70(5)$ | $0.30$
J1300+1240 | $-0.08(8)$ | $0.35$
J1400$-$1431 | $2.90(12)$ | $0.48$
J1544+4937 | $-0.6(5)$ | $0.046$
J1552+5437 | $-2.80(16)$ | $0.50$
J1640+2224 | $0.42(4)$ | $0.038$
J1658+3630 | $1.7(2)$ | $0.44$
J1713+0747 | $-0.58(15)$ | $0.052$
J1730$-$2304 | $3.5(5)$ | $1.9$
J1738+0333 | $-6.3(6)$ | $4.3$
J1744$-$1134 | $0.09(10)$ | $0.046$
J1853+1303 | $1.0(3)$ | $0.34$
J1857+0943 | $-1.7(6)$ | $1.5$
J1911$-$1114 | $-0.1(6)$ | $2.2$
J1918$-$0642 | $-1.6(14)$ | $2.7$
J1923+2515 | $2.00(11)$ | $0.33$
J1944+0907 | $-9.0(5)$ | $5.4$
J1955+2908 | -7(2) | $23$
J2043+1711 | $0.90(9)$ | $0.23$
J2051$-$0827 | $-7.1(3)$ | $2.9$
J2145$-$0750 | $0.52(4)$ | $0.069$
J2222$-$0137 | $0.1(2)$ | $0.12$
J2302+4442 | $-9.0(4)$ | $6.2$
J2317+1439 | $-0.19(15)$ | $1.3$

###  4.3 Structure functions

The IISM turbulence is known to often follow a Kolmogorov spectrum (Armstrong
et al. 1995). To characterise the turbulence in our dataset, we investigated
the structure function of the DM time series, $D_{\text{DM}}$, which gives
rise to a statistical correlation between pairs of DMs at particular time lags
$\tau$. It is defined as follows (see, e.g. You et al. 2007):

$D_{\text{DM}}(\tau)=\left<\left[\text{DM}(t+\tau)-\text{DM}(t)\right]^{2}\right>,$
(4)

where the angle brackets indicate the ensemble average. For Kolmogorov
turbulence, the structure function of the DM is expected to be a power law of
the form:

$D_{\text{DM}}(\tau)\propto\tau^{5/3}.$ (5)

In our analysis, we evaluated the structure function at equal intervals on a
logarithmic scale. Figures 5 and 6 show the structure functions for the DM
time series presented in Figs. 2 and 3, that is, all pulsars that are not
significantly affected by the Solar wind. The amplitude of the Kolmogorov
structure function fits (see below) at a time lag of 1000 days is given in
Table 2.

To investigate whether our structure functions are compatible with the
prediction arising from the assumption of a Kolmogorov spectrum, we simulated
1000 DM time series per pulsar with the same sampling as our observations
assuming a Kolmogorov power spectrum (i.e. with a spectral index of -8/3) and
calculated the structure functions. From this set of simulated structure
functions we took the median value at each time lag as well as the 68%
confidence interval.

To get an amplitude for the simulation, we first ran a least-squares routine
to fit the amplitude of a Kolmogorov power law to the simulated structure
function, weighting each sample by its inverse variance. Then we fit a
Kolmogorov power law plus a constant (to account for the white noise) to the
observed structure function, adding the relative uncertainties of the data and
the simulations in quadrature to obtain a relative uncertainty to use in the
fit. This choice of weights in the fit accounts for the high red-noise
uncertainty at large lags, which is accounted for in our simulations but not
in the data-derived structure function. Using the quotient of the two fit
amplitudes, we re-scaled the simulated structure function to match the
observed one. We then subtracted the white-noise level from the observed
structure function, such that the resulting structure function only contained
the effects of turbulent processes.

In Figure 7, we present structure functions of the DM time series that are
affected by the Solar wind (Fig. 1). As we aim to quantify the IISM
turbulence, we mitigated the effect of the Solar wind by sampling the DM time
series when the pulsar was furthest from the Sun. We then interpolated between
those samples using a cubic spline. The structure function was only calculated
for time lags greater than half a year as the short-term variations are
underestimated from the interpolation.

From our simulations we conclude that all our structure functions are
consistent with a Kolmogorov turbulence spectrum. Many sources show a decrease
(e.g. PSRs J1300+1240 and J1125+7819) or dips (e.g. PSRs J1022+1001 and
J2145$-$0750) in the structure function at large lags. The reason behind this
is that the large lags strongly depend on the start and end date of the
observations, as those define a very specific realisation of the stochastic
red-noise process. This is also reflected in the growing uncertainties of our
simulations at large lags and can be seen in Keith et al. (2013, PSRs
J0437$-$4715 or J0711$-$6830). For sources with clear linear trends in DM,
like PSRs J0740+6620 or J1400$-$1431, the structure function tends to look
steeper than Kolmogorov turbulence, while still being consistent with it. This
is expected as the structure function of a linear gradient in DM has a
spectral index of 2 (compared to a spectral index of 5/3 for Kolmogorov
turbulence). As discussed in Sect. 4.2, such a linear gradient can be
explained by the pulsars’ motion along the line of sight (Lam et al. 2016).
PSR J2317+1439 shows an abrupt reversal in the slope of the DM time series
(see Fig. 1), which is unusual for a Kolmogorov process. Further monitoring of
this source will reveal whether the structure function will remain consistent
with Kolmogorov turbulence for time lags above $\sim$1000 days.

Figure 5: Structure functions for the DM time series presented in Fig. 2. The
dashed lines indicate the median and 68% confidence interval of the Kolmogorov
structure function simulations (see text for details). The solid black
vertical line indicates $\tau=1\,\text{yr}$. The white noise level has been
subtracted from the data and the Kolmogorov model. Arrows indicate upper
limits for measurements consistent with zero at the 1-$\sigma$ level.  Figure
6: Same as Fig. 5, but for the DM time series presented in Fig. 3.  Figure 7:
Structure functions for the DM time series presented in Fig. 1 after
mitigation of the Solar wind signal. The representation is the same as in Fig.
5.

###  4.4 Comparison to PTA results

As the mitigation of the IISM is an important consideration in PTA
gravitational-wave analysis, DM time series of PTA pulsars are studied
extensively (e.g. Keith et al. 2013; Desvignes et al. 2016; Jones et al.
2017), although these use much lower cadence and are obtained at much higher
frequencies (typically $\sim$300-2000 MHz). The published PTA DM time series
present datasets that mostly end around the start of our dataset, making a
direct comparison of the DM time series difficult.

Keith et al. (2013) used observations from three frequency bands centred at
$\sim$700 MHz, $\sim$1400 MHz and $\sim$3100 MHz. From these data, they
derived one DM estimate every three months for the 20 pulsars in their sample.
As this work used the Parkes Pulsar Timing Array (PPTA) observations from the
Parkes telescope (which is located in the southern hemisphere), there are only
seven pulsars from their study that we also observed – all of which are at a
low declination where the sensitivity of LOFAR decreases dramatically (e.g.
Noutsos et al. 2015). Keith et al. (2013) do not see clear deviations from the
simple symmetric Solar wind model they applied, but this is expected as their
DM precision is not high enough. However, for some pulsars (e.g. PSR
J1857+0943), their DM precision is comparable to the one achieved with LOFAR.
Keith et al. (2013) also give the value of the structure function at a time
lag of 1000 days, which is consistent with our results within the red-noise
uncertainty.

Desvignes et al. (2016) present timing results for 42 MSPs observed with
telescopes of the European Pulsar Timing Array (EPTA). The observations were
taken across a range of frequencies between 350 MHz and 2.6 GHz, with the
majority being taken around 1.4 GHz. Of the pulsars in their study, 20 are in
common with this work. In their paper, they do not present DM time series, but
instead used the temponest software to model the DM variations as a second-
order polynomial plus a spectral noise model. This makes a comparison of
individual measurements difficult, but the large-scale trends can be compared.
Whether or not the LOFAR data provide more information on the DM variability
is again strongly pulsar dependent. For example, we detect a very clear DM
signal in PSR J2317+1439 (see Fig. 1) while Desvignes et al. (2016) only find
a DM trend at the $1\sigma$-level, consistent with the overall trend we see in
the early part of our dataset.

Jones et al. (2017) present data from the North American Nanohertz Observatory
for Gravitational Waves (NANOGrav) collaboration, analysing the DM time
evolution of 37 MSPs at frequencies between 300 MHz and 2.4 GHz, 18 of which
we analysed in this study. Again, the question of which instrument provides
greater precision in DM measurements is pulsar dependent: LOFAR provides a
much higher precision for PSR J1012+5307, whereas the higher-frequency
NANOGrav data provide a much higher precision for PSR J1857+0943 (PSR B1855+09
in Jones et al. 2017).

Overall, we find no inconsistency between our results and the highlighted
recent PTA publications. However, this is partly caused by the fact that the
datasets are difficult to compare due to a lack of overlap in observing
epochs, and differences in data representation.

###  4.5 Consequences for PTAs

In pulsar timing experiments, the actual ToA is usually not the measure of
interest. Instead, the DM-corrected, infinite-frequency ToA $T_{\infty}$ is
the relevant measure when trying to measure IISM-independent effects. As
variations in the DM at a relevant magnitude are very common, they usually
have to be accounted for (Verbiest et al. 2016). In the following, we will
discuss different approaches to apply corrections of the dispersive delays to
ToAs at 1.4 GHz, with a timing precision of
$\sigma_{\text{ToA}}=\text{$1\text{\,}\mathrm{\SIUnitSymbolMicro s}$}$.

There are two basic approaches to consider. Approach I would be to calculate
the in-band DM (as in this paper) and apply the according dispersive delay
(using Eq. 1) to the timing ToAs. The uncertainty on $T_{\infty}$ would then
be

$\sigma_{\infty}=\sqrt{\sigma_{\text{ToA}}^{2}+\sigma_{\text{disp}}^{2}}.$ (6)

To avoid the dispersive correction of being the dominant factor in this
example ($\sigma_{\text{ToA}}=\text{$1\text{\,}\mathrm{\SIUnitSymbolMicro
s}$}$ at 1.4 GHz), the uncertainty on the DM would have to be lower than
$5\times 10^{-4}$ cm-3 pc, which would not be achievable from in-band
measurements at 1.4 GHz with a bandwidth of $\sim$300 MHz, but potentially
from simultaneous lower-frequency DM measurements as presented in this paper.
Additionally, high-cadence DM time series could be smoothed to increase the
precision of the dispersive corrections if the DM variations are smooth.

Approach II would be to use ToAs from multiple frequency bands and fit for the
DM and $T_{\infty}$ simultaneously. This has the advantage of using the entire
available bandwidth, but the disadvantage that exactly simultaneous
observations may not be available, in which case short-term signatures from
any timing parameter (e.g. parameters describing the binary motion) have to be
ruled out. Also, non-simultaneous observations can have a significantly
different DM if they are too far separated in time (see Lam et al. 2015).

Lee et al. (2014) investigated how the choice of the observing frequencies
$\nu_{i}$ and the corresponding ToA uncertainties $\sigma_{i}$ affect the
precision of the infinite-frequency ToAs. Equation 12 of their paper shows
this relation in the case of only two observing frequencies:

$\left<\delta
T_{\infty}^{2}\right>=\frac{\sigma_{1}^{2}\nu_{1}^{4}+\sigma_{2}^{2}\nu_{2}^{4}}{(\nu_{1}^{2}-\nu_{2}^{2})^{2}}.$
(7)

The $\sigma_{\infty}$ can be computed as the square root of the above
expression. While the special case of only two ToAs is simplified, the general
conclusions also hold for cases with additional observing frequencies. Most
importantly, the frequency band with the largest $Q=\sigma^{2}\nu^{4}$
dominates the uncertainty in $T_{\infty}$. Due to the strong frequency
dependence of $Q$, the highest frequency has to be the most-precisely measured
in order to benefit from the additional bandwidth. If $\nu_{1}$ is our most
precise timing frequency and observations at $\nu_{1}$ and $\nu_{2}$ are used
to calculate $T_{\infty}$ ($\sigma_{2}>\sigma_{1}$), Fig. 8 illustrates that
$\nu_{2}$ should be smaller than $\nu_{1}$ to achieve a timing precision close
to $\sigma_{1}$. If $\nu_{2}$ is the higher frequency, $\sigma_{\infty}$ is
limited by $\sigma_{2}$. In PTA data, the smallest ToA uncertainty is often
obtained at 1.4 GHz, so observations at lower frequencies should be used to
correct for DM variability. In the case of the LOFAR observations presented in
this paper, a frequency-integrated ToA uncertainty of
$\sim$$100\text{\,}\mathrm{\SIUnitSymbolMicro s}$ or better would be required
to not dominate DM corrections at the $1\text{\,}\mathrm{\SIUnitSymbolMicro
s}$-level. This condition is fulfilled for all pulsars in our sample except
for PSR J0621+1001 (see Table 1), which is a pulsar that is usually timed at
an RMS much worse than $1\text{\,}\mathrm{\SIUnitSymbolMicro s}$ (see e.g.
Desvignes et al. 2016). This shows that Approach II should be preferred over
Approach I, if it is applicable.

Figure 8: Illustration of Eq. 7 (Eq. 12 in Lee et al. 2014) for an observation
at two frequency bands. The uncertainty of the infinite-frequency ToA is
plotted as a function of $\nu_{2}$, with $\nu_{2}$ providing less precise ToAs
than $\nu_{1}$ ($\sigma_{2}=5.0\sigma_{1}$). All quantities are given relative
to the frequency and ToA uncertainty of Band 1. On the left side of the plot
($\nu_{2}\ll\nu_{1}$), $\sigma_{\infty}$ approaches $\sigma_{1}$, whereas on
the right side ($\nu_{2}\gg\nu_{1}$) it approaches $\sigma_{2}$. This shows
that the precision of the higher-frequency band is the limiting factor on the
precision of the infinite-frequency ToA.

To illustrate the superiority of Approach II over Approach I, we will go
through an example of correcting for the dispersive delay in 1.4 GHz
observations using LOFAR. We assume $\nu_{1}=1.4\,\text{GHz}$,
$\sigma_{1}=\text{$1\text{\,}\mathrm{\SIUnitSymbolMicro s}$}$,
$\nu_{2}=150\,\text{MHz}$ and
$\sigma_{2}=\text{$100\text{\,}\mathrm{\SIUnitSymbolMicro s}$}$. We note that
a ToA precision of $100\text{\,}\mathrm{\SIUnitSymbolMicro s}$ is a worst-case
scenario. Using Approach II, we get
$\sigma_{\infty}=\text{$1.5\text{\,}\mathrm{\SIUnitSymbolMicro s}$}$. For
Approach I, we assume ten ToAs around $\nu_{2}$ across a bandwidth of 75 MHz,
each with an uncertainty of
$\sqrt{10}\times\text{$100\text{\,}\mathrm{\SIUnitSymbolMicro s}$}$ and use
Eq. 10 of Lee et al. (2014) to compute the DM uncertainty. The resulting
$\sigma_{\infty}$ is $3.8\text{\,}\mathrm{\SIUnitSymbolMicro s}$, which is
worse by a factor of 2.5.

A major caveat of using LOFAR observations for DM corrections at higher
frequencies lies in the spectra of pulsars: Pulsars with steep spectra are
bright at low frequencies and allow for highly precise DM corrections, but are
fainter at high frequencies, such that the timing precision may be so low that
DM corrections are less relevant. On the other hand, pulsars with flatter
spectra tend to have a high ToA precision at high frequencies, but are very
faint at low frequencies due to the steep spectral index of the sky background
noise, which is of order 2.5 (Lawson et al. 1987). For example, the pulsar
with the best DM precision in our sample is PSR J0034$-$0534, but its timing
RMS at 1.4 GHz is only average ($4.27\text{\,}\mathrm{\SIUnitSymbolMicro s}$,
see Perera et al. 2019). The well-timed PSR J1713+0747, however, is so faint
we cannot even detect it in GLOW data and achieve a DM precision an order of
magnitude worse than for PSR J0034$-$0534.

It therefore follows that Approach I using LOFAR data would currently not
improve the PTA timing precision for most sources, which could be solved by
increasing the timing precision of steep-spectrum sources at 1.4 GHz or
improving the DM precision of flat-spectrum sources at low frequencies. The
latter could be achieved by greatly increasing the integration time with the
LOFAR Core or using more sensitive telescopes like the upcoming Square
Kilometre Array phase-one low-frequency bands (SKA1-LOW; e.g. Janssen et al.
2015). Using a slightly higher frequency might also help if the spectrum of
the pulsar is very flat – the optimal frequency for the DM measurement is
strongly pulsar dependent. As the timing precision of steep-spectrum PTA
sources continues to increase due to improved telescopes, the significance of
highly precise DM time series such as those presented in this work will become
crucial for DM corrections. However, Approach II should already now improve
the timing precision of PTAs for many sources, as it provides more precise
corrections for the dispersive delays. We reserve such an analysis for a
future paper, however.

###  4.6 Data access

Our ToAs, timing models, templates and DM time series are available online on
Zenodo888https://zenodo.org/deposit/4290012
DOI: 10.5281/zenodo.4290012. The raw observations taken with the LOFAR Core
can be accessed via the LOFAR Long Term Archive999https://lta.lofar.eu/. The
raw GLOW data are available upon request.

##  5 Conclusions

We present low-frequency DM time series for 36 MSPs over up to 7.1 years,
obtained from observations with the LOFAR Core and the individual GLOW
telescopes. Except for the pulsars with very high DM uncertainty (i.e. greater
than $\sim 10^{-4}$ cm-3 pc), all pulsars show significant variations in DM.
Twelve pulsars show a clear Solar wind signal in their DM time series that can
usually not be modelled by a spherically symmetric electron content with
constant amplitude, which is often used in pulsar timing. All of the IISM-
related DM variations we present are consistent with a Kolmogorov turbulence
spectrum.

Acknowledging the caveat that LOFAR often provides a high DM precision for
pulsars that are poorly timed at higher frequencies, and vice versa, we show
that our LOFAR DM monitoring could be used to correct variations of the
dispersive delays in higher-frequency observations from PTAs. We do not find
evidence for a frequency-dependent DM, so we expect the impact of this effect
to be limited.

###### Acknowledgements.

We like to thank W. A. Coles and M. Lam for useful discussions. JPWV
acknowledges support by the Deutsche Forschungsgemeinschaft (DFG) through the
Heisenberg programme (Project No. 433075039). CT is a Veni fellow (project
number 016.Veni.192.086), partly financed by the Dutch Re-search Council
(NWO). RPB acknowledges support of the European Research Council, under the
European Union’s Horizon 2020 research and innovation program (grant agreement
No. 715051; Spiders). JvL acknowledges funding from the European Research
Council under the European Union’s Seventh Framework Programme (FP/2007-2013)
/ ERC Grant Agreement n. 617199 (‘ALERT’), and from Vici research programme
‘ARGO’ with project number 639.043.815, financed by the Dutch Research Council
(NWO). Part of this work is based on data obtained with the International
LOFAR Telescope under project codes LC0_011, LC1_027, LC2_010, LT3_001,
LC4_004, LT5_003, LC9_041 and LT10_004. LOFAR (van Haarlem et al. 2013) is the
Low Frequency Array designed and constructed by ASTRON. It has observing, data
processing, and data storage facilities in several countries, that are owned
by various parties (each with their own funding sources), and that are
collectively operated by the ILT foundation under a joint scientific policy.
The ILT resources have benefitted from the following recent major funding
sources: CNRS-INSU, Observatoire de Paris and Université d’Orléans, France;
BMBF, MIWF-NRW, MPG, Germany; Science Foundation Ireland (SFI), Department of
Business, Enterprise and Innovation (DBEI), Ireland; NWO, The Netherlands; The
Science and Technology Facilities Council, UK. This paper uses data obtained
with the German LOFAR stations, during station-owners time and ILT time
allocated under project codes LC0_014, LC1_048, LC2_011, LC3_029, LC4_025,
LT5_001, LC9_039 and LT10_014. We made use of data from the Effelsberg (DE601)
LOFAR station funded by the Max-Planck-Gesellschaft; the Unterweilenbach
(DE602) LOFAR station funded by the Max-Planck-Institut für Astrophysik,
Garching; the Tautenburg (DE603) LOFAR station funded by the State of
Thuringia, supported by the European Union (EFRE) and the Federal Ministry of
Education and Research (BMBF) Verbundforschung project D-LOFAR I (grant
05A08ST1); the Potsdam (DE604) LOFAR station funded by the Leibniz-Institut
für Astrophysik, Potsdam; the Jülich (DE605) LOFAR station supported by the
BMBF Verbundforschung project DLOFAR I (grant 05A08LJ1); and the Norderstedt
(DE609) LOFAR station funded by the BMBF Verbundforschung project D-LOFAR II
(grant 05A11LJ1). The observations of the German LOFAR stations were carried
out in stand-alone GLOW mode, which is technically operated and supported by
the Max-Planck-Institut für Radioastronomie, the Forschungszentrum Jülich and
Bielefeld University. We acknowledge support and operation of the GLOW
network, computing and storage facilities by the FZ-Jülich, the MPIfR and
Bielefeld University and financial support from BMBF D-LOFAR III (grant
05A14PBA) and D-LOFAR IV (grant 05A17PBA), and by the states of Nordrhein-
Westfalia and Hamburg. We acknowledge the work of A. Horneffer in setting up
the GLOW network and initial recording machines.

## References

  * Armstrong et al. (1995) Armstrong, J. W., Rickett, B. J., & Spangler, S. R. 1995, ApJ, 443, 209
  * Arzoumanian et al. (2018) Arzoumanian, Z., Baker, P. T., Brazier, A., et al. 2018, ApJ, 859, 47
  * Arzoumanian et al. (2015) Arzoumanian, Z., Brazier, A., Burke-Spolaor, S., et al. 2015, ApJ, 813, 65
  * Burke-Spolaor et al. (2019) Burke-Spolaor, S., Taylor, S. R., Charisi, M., et al. 2019, A&A Rev., 27, 5
  * Cordes et al. (2016) Cordes, J. M., Shannon, R. M., & Stinebring, D. R. 2016, ApJ, 817, 16
  * Demorest (2011) Demorest, P. B. 2011, MNRAS, 416, 2821
  * Demorest et al. (2013) Demorest, P. B., Ferdman, R. D., Gonzalez, M. E., et al. 2013, ApJ, 762, 94
  * Desvignes et al. (2016) Desvignes, G., Caballero, R. N., Lentati, L., et al. 2016, MNRAS, 458, 3341
  * Donner et al. (2019) Donner, J. Y., Verbiest, J. P. W., Tiburzi, C., et al. 2019, A&A, 624, A22
  * Hellings & Downs (1983) Hellings, R. W. & Downs, G. S. 1983, ApJ, 265, L39
  * Hobbs & Dai (2017) Hobbs, G. & Dai, S. 2017, National Science Review, 4, 707
  * Hobbs et al. (2010) Hobbs, G., Lyne, A. G., & Kramer, M. 2010, MNRAS, 402, 1027
  * Hobbs et al. (2004) Hobbs, G., Lyne, A. G., Kramer, M., Martin, C. E., & Jordan, C. 2004, MNRAS, 353, 1311
  * Hobbs et al. (2006) Hobbs, G. B., Edwards, R. T., & Manchester, R. N. 2006, MNRAS, 369, 655
  * Hotan et al. (2005) Hotan, A. W., Bailes, M., & Ord, S. M. 2005, MNRAS, 362, 1267
  * Hotan et al. (2004) Hotan, A. W., van Straten, W., & Manchester, R. N. 2004, PASA, 21, 302
  * Jammalamadaka & SenGupta (2001) Jammalamadaka, S. R. & SenGupta, A. 2001, Topics in Circular Statistics (New Jersey: World Scientific)
  * Janssen et al. (2015) Janssen, G., Hobbs, G., McLaughlin, M., et al. 2015, Advancing Astrophysics with the Square Kilometre Array (AASKA14), 37
  * Jones et al. (2017) Jones, M. L., McLaughlin, M. A., Lam, M. T., et al. 2017, ApJ, 841, 125
  * Keith et al. (2013) Keith, M. J., Coles, W., Shannon, R. M., et al. 2013, MNRAS, 429, 2161
  * Kondratiev et al. (2016) Kondratiev, V. I., Verbiest, J. P. W., Hessels, J. W. T., et al. 2016, A&A, 585, A128
  * Kulkarni (2020) Kulkarni, S. R. 2020, arXiv e-prints, arXiv:2007.02886
  * Künkel (2017) Künkel, L. 2017, Master’s thesis, Bielefeld University, Faculty of Physics
  * Lam et al. (2016) Lam, M. T., Cordes, J. M., Chatterjee, S., et al. 2016, ApJ, 819, 155
  * Lam et al. (2015) Lam, M. T., Cordes, J. M., Chatterjee, S., & Dolch, T. 2015, ApJ, 801, 130
  * Lam et al. (2018) Lam, M. T., Ellis, J. A., Grillo, G., et al. 2018, The Astrophysical Journal, 861, 132
  * Lawson et al. (1987) Lawson, K. D., Mayer, C. J., Osborne, J. L., & Parkinson, M. L. 1987, MNRAS, 225, 307
  * Lazarus et al. (2016) Lazarus, P., Karuppusamy, R., Graikou, E., et al. 2016, MNRAS, 458, 868
  * Lee et al. (2014) Lee, K. J., Bassa, C. G., Janssen, G. H., et al. 2014, MNRAS, 441, 2831
  * Levin et al. (2016) Levin, L., McLaughlin, M. A., Jones, G., et al. 2016, ApJ, 818, 166
  * Lorimer & Kramer (2005) Lorimer, D. R. & Kramer, M. 2005, Handbook of Pulsar Astronomy (Cambridge University Press)
  * Manchester et al. (2005) Manchester, R. N., Hobbs, G. B., Teoh, A., & Hobbs, M. 2005, AJ, 129, 1993
  * Maron et al. (2000) Maron, O., Kijak, J., Kramer, M., & Wielebinski, R. 2000, A&AS, 147, 195
  * Noutsos et al. (2015) Noutsos, A., Sobey, C., Kondratiev, V. I., et al. 2015, A&A, 576, A62
  * Perera et al. (2019) Perera, B. B. P., DeCesar, M. E., Demorest, P. B., et al. 2019, MNRAS, 490, 4666
  * Sanidas et al. (2019) Sanidas, S., Cooper, S., Bassa, C. G., et al. 2019, A&A, 626, A104
  * Stappers et al. (2011) Stappers, B. W., Hessels, J. W. T., Alexov, A., et al. 2011, A&A, 530, A80
  * Tauris et al. (2015) Tauris, T. M., Kaspi, V. M., Breton, R. P., et al. 2015, in Advancing Astrophysics with the Square Kilometre Array (AASKA14), 39
  * Taylor (1992) Taylor, J. H. 1992, Philos. Trans. Roy. Soc. London A, 341, 117
  * Tiburzi (2018) Tiburzi, C. 2018, PASA, 35, e013
  * Tiburzi et al. (2019) Tiburzi, C., Verbiest, J. P. W., Shaifullah, G. M., et al. 2019, MNRAS, 487, 394
  * van Haarlem et al. (2013) van Haarlem, M. P., Wise, M. W., Gunst, A. W., et al. 2013, A&A, 556, A2
  * van Straten & Bailes (2011) van Straten, W. & Bailes, M. 2011, PASA, 28, 1
  * van Straten et al. (2012) van Straten, W., Demorest, P., & Oslowski, S. 2012, Astronomical Research and Technology, 9, 237
  * Verbiest et al. (2009) Verbiest, J. P. W., Bailes, M., Coles, W. A., et al. 2009, MNRAS, 400, 951
  * Verbiest et al. (2016) Verbiest, J. P. W., Lentati, L., Hobbs, G., et al. 2016, MNRAS, 458, 1267
  * You et al. (2007) You, X.-P., Hobbs, G., Coles, W., et al. 2007, MNRAS, 378, 493

  *[LOFAR]: LOw Frequency ARray
  *[MPIfR]: Max-Planck Institut für Radioastronomie
  *[DM]: dispersion measure
  *[PTAs]: pulsar timing array
  *[ToAs]: time of arrival
  *[ToA]: time of arrival
  *[MSPs]: millisecond pulsar
  *[PTA]: pulsar timing array
  *[IISM]: ionised interstellar medium
  *[S/N]: signal-to-noise ratio
  *[DMs]: dispersion measure
  *[]:
  *[ILT]: International   () Telescope
  *[PulP]: the pulsar pipeline
  *[GLOW]: German LOng Wavelength
  *[LuMP]:  und   () Pulsare
  *[EPTA]: European Pulsar Timing Array
  *[PPTA]: Parkes Pulsar Timing Array
  *[NANOGrav]: North American Nanohertz Observatory for Gravitational Waves
  *[IPTA]: International Pulsar Timing Array
  *[RFI]: radio frequency interference


# F-PABEE: Flexible-patience-based Early Exiting for Single-label and Multi-
label text Classification Tasks

###### Abstract

Computational complexity and overthinking problems have become the bottlenecks
for pre-training language models (PLMs) with millions or even trillions of
parameters. A Flexible-Patience-Based Early Exiting method (F-PABEE) has been
proposed to alleviate the problems mentioned above for single-label
classification (SLC) and multi-label classification (MLC) tasks. F-PABEE makes
predictions at the classifier and will exit early if predicted distributions
of cross-layer are consecutively similar. It is more flexible than the
previous state-of-the-art (SOTA) early exiting method PABEE because it can
simultaneously adjust the similarity score thresholds and the patience
parameters. Extensive experiments show that: (1) F-PABEE makes a better
speedup-accuracy balance than existing early exiting strategies on both SLC
and MLC tasks. (2) F-PABEE achieves faster inference and better performances
on different PLMs such as BERT and ALBERT. (3) F-PABEE-JSKD performs best for
F-PABEE with different similarity measures.

Index Terms—  F-PABEE, PABEE, Early Exiting, Multi-label Classification,
Single-label Classification

##  1 Introduction

Fine-tuning PLMs has become the de-facto paradigm in natural language
processing [1], due to the amazing performance gains on a wide range of
natural language processing tasks [2, 3, 4, 5, 6]. Despite SOTA performances,
BERT [7] and its variants [8, 9, 10, 11] still face significant application
challenges: cumbersome computation and overthinking problems due to huge
parameters and deep models. Early exiting attracts much attention as an input-
adaptive method to speed up inference [12]. Early exiting installs a
classifier at each transformer layer to evaluate the predictions and will exit
when meeting the criterion. Three different early exiting strategies exist:
(1) The confidence-based strategy evaluates the predictions based on specific
confidence measurements. (2) The learned-based strategy learns a criterion for
early exiting. (3) The patience-based strategy exits when consecutive
classifiers make the exact predictions. Among them, the patience-based
strategy PABEE [13] achieves SOTA results.

We raise two issues for the current SOTA strategy: (1) PABEE faces a
limitation for application: it can not flexibly adjust the speedup ratio on a
given task and fixed patience parameter, mainly caused by a strict cross-layer
comparison strategy. Thus, we wonder whether we can combine PABEE with a
softer cross-layer comparison strategy. (2) Current early exiting strategies
mainly focus on SLC tasks, while the MLC tasks are neglected. So can they
speed up MLC tasks?

Therefore, we propose a Flexible-Patience-Based Early Exiting method (F-PABEE)
to address the above issues. F-PABEE makes predictions at each classifier and
will exit early if the current layer and the last few layers have similar
(similarity score less than a threshold) predicted distributions. F-PABEE can
be seen as a natural extension of PABEE and is more flexible since it can
achieve better speed-accuracy tradeoffs by adjusting the similarity score
thresholds and patience parameters. It can also extend to MLC tasks
effortlessly.

Our contributions are summarized as follows: (1) We propose F-PABEE, a novel
and effective inference mechanism that is flexible in adjusting the speedup
ratios of PLMs. (2) The results show that our method can accelerate inference
effectively while maintaining good performances across different SLC and MLC
tasks. (3) We are the first to investigate the early exiting of MLC tasks, and
F-PABEE is suitable for this type of task.

Fig. 1: Inference procedure of PABEE and F-PABEE, $C_{i}$ is the classifier,
$thre$ is threshold, $P_{0}$ is pre-defined patience.

##  2 Related works

###  2.1 Static inference approach

The static inference approach compresses the heavy model into a smaller one,
including pruning, knowledge distillation, quantization, and weight sharing
[14, 15, 16]. For example, HeadPrune [17] ranks the attention heads and prunes
them to reduce inference latency. PKD [18] investigates the best practices of
distilling knowledge from BERT into smaller-sized models. I-BERT [19] performs
an end-to-end BERT inference without any floating point calculation. ALBERT
[8] shares the cross-layer parameters. [20], [21] and [22] distills knowledge
from the larger BERT teacher model for improving the performances of student
networks which are learned with neural architecture search. Note that the
static models are still in the form of deep neural networks with multiple
stacked layers. The computational path is invariable for all examples in the
inference process, which is not flexible.

###  2.2 Dynamic early exiting

Orthogonal to the static inference approach, early exiting dynamically adjusts
hyper-parameters in response to changes in request traffic. It does not need
to make significant changes to the original model structure or weight bits,
nor does it need to train different teacher-student learning networks, which
saves computing resources [23].

There are mainly three groups of dynamic early exiting strategies. The first
type is confidence-based early exiting [24]. For example, BranchyNet [25],
FastBERT [26], and DeeBERT [27] calculate the entropy of the prediction
probability distribution to estimate the confidence of classifiers to enable
dynamic early exiting. Shallow-deep [28] and RightTool [29] leverage the
maximum of the predicted distribution as the exiting signal. The second type
is the learned-based exiting, such as BERxiT [30] and CAT [31]. They learn a
criterion for early exiting. The third type is patience-based early exiting,
such as PABEE [13], which stops inference and exits early if the classifiers’
predictions remain unchanged for pre-defined times. Among them, patience-based
PABEE achieves SOTA performance. However, PABEE suffers from too strict cross-
layer comparison, and the applications on MLC tasks are neglected. There are
also literature focusing on improving the training of multi-exit BERT, like
LeeBERT [32] and GAML-BERT [33].

F-PABEE is a more flexible extension to PABEE, which can simultaneously adjust
the confidence thresholds and patience parameters to meet different
requirements. In addition, it outperforms other existing early exiting
strategies on both SLC and MLC tasks.

###  2.3 Training of multi-exit backbones

The literature on early exiting focuses more on the design of early exiting
strategies, thus neglect the advances of multi-exit backbones’ training
methods. LeeBERT [32] employs an adaptive learning method for training
multiple exits. GAML-BERT [33] enhances the training of multi-exit backbones
by a mutual learning approach.

| CoLA | MNLI | MRPC | QNLI | QQP | RTE | SST-2
---|---|---|---|---|---|---|---
| score | speedup | score | speedup | score | speedup | score | speedup | score | speedup | score | speedup | score | speedup
BERT base | 54.2 | 0% | 83.1 | 0% | 86.8 | 0% | 89.8 | 0% | 89.2 | 0% | 69.1 | 0% | 91.3 | 0%
Fixed-Exit-3L | 0.0 | 75% | 70.0 | 75% | 75.8 | 75% | 77.4 | 75% | 81.8 | 75% | 54.7 | 75% | 81.0 | 75%
Fixed-Exit-6L | 0.0 | 50% | 79.6 | 50% | 84.7 | 50% | 85.3 | 50% | 89.3 | 50% | 68.1 | 50% | 88.6 | 50%
BranchyNet | 0.0 | 74% | 63.8 | 76% | 75.7 | 76% | 74.2 | 80% | 71.6 | 80% | 54.7 | 76% | 79.9 | 76%
0.0 | 51% | 78.3 | 53% | 83.0 | 52% | 87.1 | 47% | 89.3 | 50% | 67.4 | 47% | 88.3 | 49%
Shallow-Deep | 0.0 | 75% | 64.1 | 77% | 75.6 | 76% | 74.3 | 78% | 71.4 | 79% | 54.7 | 76% | 79.5 | 77%
0.0 | 52% | 78.2 | 51% | 82.8 | 51% | 87.2 | 49% | 89.6 | 51% | 67.2 | 48% | 88.4 | 48%
BERxiT | 0.0 | 76% | 63.5 | 76% | 75.6 | 76% | 73.3 | 78% | 68.2 | 80% | 55.3 | 77% | 79.5 | 76%
12.3 | 52% | 78.4 | 51% | 82.9 | 51% | 87.0 | 48% | 89.1 | 49% | 67.3 | 47% | 88.3 | 49%
PABEE | 0.0 | 75% | 63.9 | 77% | 75.8 | 75% | 73.6 | 81% | 68.6 | 82% | 55.8 | 75% | 79.9 | 77%
0.0 | 50% | 78.9 | 52% | 83.1 | 53% | 87.2 | 46% | 89.6 | 49% | 67.7 | 46% | 88.7 | 48%
F-PABEE | 0.0 | 75% | 66.9 | 72% | 81.5 | 77% | 76.2 | 75% | 79.6 | 82% | 56.0 | 76% | 80.5 | 76%
13.6 | 52% | 83.9 | 53% | 87.3 | 53% | 88.6 | 54% | 90.8 | 49% | 68.1 | 47% | 92.3 | 48%

Table 1: Experimental results of different early exiting methods with BERT
backbone on the GLUE benchmark.

##  3 Flexible patience-based early exiting

###  3.1 Inference procedure for SLC and MLC tasks

The inference procedure of F-PABEE is shown in Fig 1(b), which is an improved
version of PABEE (Fig 1(a)), where $L_{i}$ is the transformer block of the
model, $n$ is the number of transformer layers, $C_{i}$ is the inserted
classifier layer, $s$ is the cross-layer similarity score, $thre$ is the
similarity score threshold, $P_{0}$ is the pre-defined patience value in the
model.

The input sentences are first embedded as the vector:

$h_{0}=\text{Embedding}(x).$ (1)

The vector is then passed through transformer layers ($L_{1}...L_{n}$) to
extract features and compute its hidden state $h$. After which, we use
internal classifiers ($C_{1}...C_{n}$), which are connected to each
transformer layer to predict probability $p$:

$p_{i}=C_{i}(h_{i})=C_{i}(L_{i}(h_{i-1})).$ (2)

We denote the similarity score between the prediction results of layer $i-1$
and $i$ as $s(p_{i-1},p_{i})$ ($s(p_{i-1},p_{i})\in\mathbf{R}$). The smaller
the value of $s(p_{i-1},p_{i})$, the prediction distributions are more
consistent with each other. The premise of the model’s early exit is that the
comparison scores between successive layers are relatively small; The
similarity threshold $thre$ is a hyper-parameter. We use $pat_{i}$ to store
the times that the cross-layer comparison scores are consecutively less than
the threshold $thre$ when the model reaches current layer $i$:

$pat_{i}=\left\\{\begin{array}[]{rcl}pat_{i-1}+1&&s(p_{i-1},p_{i})<thre\\\
0&&s(p_{i-1},p_{i})>=thre\end{array}\right\\}$ (3)

If $s(p_{i-1},p_{i})$ is less than the similarity score threshold $thre$, then
increase the patience counter by 1. Otherwise, reset the patience counter to
0. This process is repeated until $pat$ reaches the pre-defined patience value
$P_{0}$. The model dynamically stops inference and exits early. However, if
this condition is never met, the model uses the final classifier layer to make
predictions. This way, the model can stop inference early without going
through all layers.

###  3.2 Similarity measures for SLC and MLC tasks

Under the framework of F-PABEE, we can adopt different similarity measures for
predicted probability distributions. This work uses the knowledge distillation
objectives as the similarity measures [34]. When the model reaches the current
layer $l$, for SLC tasks, we compare a series of similarity measures of
F-PABEE, denoted as:

F-PABEE-KD: It adopts the knowledge distillation objective from probability
mass distribution $p^{l-1}$ to $p^{l}$:

$s(p^{l-1},p^{l})=-\sum_{j=1}^{k}p_{j}^{l-1}log(p_{j}^{l});$ (4)

F-PABEE-ReKD: It adopts the knowledge distillation objective in the reverse
direction, from probability mass distribution $p^{l}$ to $p^{l-1}$:

$s(p^{l},p^{l-1})=-\sum_{j=1}^{k}p_{j}^{l}log(p_{j}^{l-1});$ (5)

F-PABEE-SymKD: It adopts a symmetrical knowledge distillation objective:

$SymKD=s(p^{l-1},p^{l})+s(p^{l},p^{l-1});$ (6)

F-PABEE-JSKD: It adopts another symmetrical distillation objective, similar to
Jenson-Shannon divergence:

$JSKD=\frac{1}{2}s(p^{l-1},\frac{p^{l-1}+p^{l}}{2})+\frac{1}{2}s(p^{l},\frac{p^{l-1}+p^{l}}{2})$
(7)

In addition, for MLC tasks, we transform them into multiple binary
classification problems and sum the similarity scores of all categories, and
the formulas are denoted as:

F-PABEE-KD:

$s(p^{l-1},p^{l})=-\sum_{j=1}^{k}\sum_{i=1}^{2}p_{ji}^{l-1}log(p_{ji}^{l});$
(8)

F-PABEE-ReKD:

$s(p^{l},p^{l-1})=-\sum_{j=1}^{k}\sum_{i=1}^{2}p_{ji}^{l}log(p_{ji}^{l-1});$
(9)

The formulations of F-PABEE-SymKD and F-PABEE-JSKD for MLC tasks are similar
to those of SLC tasks.

Fig. 2: Speed-accuracy curves of F-PABEE, PABEE and BERxiT on SLC tasks with
BERT backbone.

###  3.3 Training procedure

F-PABEE is trained on SLC and MLC tasks, while the activation and loss
functions are different. For SLC tasks, we use the softmax activation function
and cross-entropy function according to the tasks. In contrast, we use the
sigmoid activation function and binary cross-entropy function for MLC tasks.

After that, we optimize the model parameters by minimizing the overall loss
function $L$, which is the weighted average of the loss terms from all
classifiers:

$L=\sum_{j=1}^{n}jL_{j}/\sum_{j=1}^{n}j$ (10)

##  4 Experiments

###  4.1 Tasks and Baselines

We evaluate F-PABEE on GLUE benchmark [35] for SLC tasks and four datasets for
MLC tasks: MixSNLPS [36], MixATS [37], AAPD [38], and Stackoverflow [39]. we
compare F-PABEE with three groups of baselines: (1) BERT-base; (2) Static
exiting; (3) Dynamic exiting methods, including BrachcyNet [40], Shallow-Deep
[28], BERxiT [30], and PABEE. Considering the flops of inferencing one with
the whole BERT as the base, the speed-up ratio is defined as the average ratio
of reduced flops due to early exiting.

###  4.2 Experimental setting

In training process, we perform grid search over the batch size of {16, 32,
128}, and learning rate of {1e-5, 2e-5, 3e-5, 5e-5} with an AdamW optimizer
[41] . The batch size in the inference process is 1. We implement F-PABEE on
the bases of HuggingFace Transformers [42]. All experiments are conducted on
two Nvidia TITAN X 24GB GPUs.

###  4.3 Overall comparisons

Fig. 3: Speed-accuracy curves of F-PABEE, PABEE and BERxiT on MLC tasks with
BERT backbone.  Fig. 4: The distribution of executed layers of MRPC and
MixSNIPS on average at different speeds (50%, 75%).  Fig. 5: Speed-accuracy
curves of F-PABEE, PABEE and BERxiT on SLC and MLC tasks with ALBERT backbone.
Fig. 6: Speed-accuracy curves of different similarity measures on SLC and MLC
tasks with BERT backbone.

In Table 1, we compare F-PABEE with other early exiting strategies. We adjust
the hyper-parameters of F-PABEE and other baselines to ensure similar speedups
with PABEE. It shows that F-PABEE balances speedup and performance better than
baselines, especially for a large speedup ratio. Moreover, we draw the score-
speedup curves for BERxiT, PABEE, and F-PABEE. It shows that F-PABEE
outperforms the baseline models on both SLC (Fig 2) and MLC tasks(Fig 3).
Furthermore, the distribution of executed layers (Fig 4) indicates that
F-PABEE can choose the faster off-ramp and achieve a better trade-off between
accuracy and efficiency by flexibly adjusting similarity score thresholds and
patience parameters.

###  4.4 Ablation studies

Ablation on different PLMs F-PABEE is flexible and can work well with other
pre-trained models, such as ALBERT. Therefore, to show the acceleration
ability of F-PABEE with different backbones, we compare F-PABEE to other early
exiting strategies with ALBERT base as the backbone. The results in Fig 5 show
that F-PABEE outperforms other early exiting strategies under different
backbones by large margins on both SLC and MLC tasks, indicating that F-PABEE
can accelerate the inference process for numerous PLMs.

Comparisons between different similarity measures We consider F-PABEE with
different similarity measures, denoted as F-PABEE-KD, F-PABEE-ReKD, F-PABEE-
SymKD, and F-PABEE-JSKD, and the results are presented in Fig 6. F-PABEE-JSKD
performs the best on both SLC and MLC tasks. We suppose that F-PABEE-JSKD is
symmetric, and the similarity discrimination is more accurate than asymmetric
measures. Therefore, it is better at determining which samples should exit at
shallow layers and which should go through deeper layers.

##  5 Conclusions

We proposed F-PABEE, a novel and efficient early exiting method that combines
PABEE with a softer cross-layer comparison strategy. F-PABEE is more flexible
than PABEE since it can achieve different speed-performance tradeoffs by
adjusting the similarity score thresholds and patience parameters. In
addition, we investigate the acceleration ability of F-PABEE with different
backbones. Moreover, we compare the performances of F-PABEE with different
similarity measures. Extensive experiments on SLC and MLC demonstrate that:
(1) F-PABEE performs better than the previous SOTA adaptive early exiting
strategies for both SLC and MLC tasks. As far as we know, we are the first to
investigate the early exiting methods for MLC tasks. (2) F-PABEE performs well
on different PLMs such as BERT and ALBERT. (3) Ablation studies show that
F-PABEE-JSKD performs best for F-PABEE with different similarity measures.

## References

  * [1] Tianyang Lin et al.,  “A survey of transformers,”  arXiv, vol. abs/2106.04554, 2021.
  * [2] Wei Zhu,  “Mvp-bert: Redesigning vocabularies for chinese bert and multi-vocab pretraining,”  ArXiv, vol. abs/2011.08539, 2020.
  * [3] Wei Zhu, Xiaofeng Zhou, Keqiang Wang, Xun Luo, Xiepeng Li, Yuan Ni, and Guo Tong Xie,  “Panlp at mediqa 2019: Pre-trained language models, transfer learning and knowledge distillation,”  in BioNLP@ACL, 2019.
  * [4] Yuhui Zuo, Wei Zhu, and Guoyong Cai,  “Continually detection, rapidly react: Unseen rumors detection based on continual prompt-tuning,”  in International Conference on Computational Linguistics, 2022.
  * [5] Wei Zhu, Peng Wang, Xiaoling Wang, Yuan Ni, and Guo Tong Xie,  “Acf: Aligned contrastive finetuning for language and vision tasks,”  ICASSP 2023 - 2023 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), 2023.
  * [6] Zhao Guo, Yuan Ni, Keqiang Wang, Wei Zhu, and Guo Tong Xie,  “Global attention decoder for chinese spelling error correction,”  in Findings, 2021.
  * [7] Jacob Devlin et al.,  “BERT: pre-training of deep bidirectional transformers for language understanding,”  ArXiv, vol. abs/1810.04805, 2018.
  * [8] Zhenzhong Lan et al.,  “ALBERT: A lite BERT for self-supervised learning of language representations,”  arXiv, vol. abs/1909.11942, 2019.
  * [9] Zhilin Yang et al.,  “Xlnet: Generalized autoregressive pretraining for language understanding,”  arXiv, vol. abs/1906.08237, 2019.
  * [10] Yinhan Liu et al.,  “Roberta: A robustly optimized BERT pretraining approach,”  arXiv, vol. abs/1907.11692, 2019.
  * [11] Wei Zhu,  “MVP-BERT: Multi-vocab pre-training for Chinese BERT,”  in Proceedings of the 59th Annual Meeting of the Association for Computational Linguistics and the 11th International Joint Conference on Natural Language Processing: Student Research Workshop, Online, Aug. 2021, pp. 260–269, Association for Computational Linguistics.
  * [12] Canwen et al. Xu,  “A survey on dynamic neural networks for natural language processing,”  arXiv, vol. abs/2202.07101, 2022.
  * [13] Wangchunshu Zhou et al.,  “BERT loses patience: Fast and robust inference with early exit,”  arXiv, vol. abs/2006.04152, 2020.
  * [14] Canwen et al.,  “Bert-of-theseus: Compressing bert by progressive module replacing,”  EMNLP, 2020.
  * [15] Victor Sanh et al.,  “Distilbert, a distilled version of BERT: smaller, faster, cheaper and lighter,”  arXiv, vol. abs/1910.01108, 2019.
  * [16] Angela Fan et al.,  “Reducing transformer depth on demand with structured dropout,”  arXiv, vol. abs/1909.11556, 2019.
  * [17] Paul Michel et al.,  “Are sixteen heads really better than one?,”  arXiv, vol. abs/1905.10650, 2019.
  * [18] Siqi Sun et al.,  “Patient knowledge distillation for BERT model compression,”  arXiv, vol. abs/1908.09355, 2019.
  * [19] Sehoon Kim et al.,  “I-BERT: integer-only BERT quantization,”  arXiv, vol. abs/2101.01321, 2021.
  * [20] Wei Zhu,  “Autonlu: Architecture search for sentence and cross-sentence attention modeling with re-designed search space,”  in Natural Language Processing and Chinese Computing, Lu Wang, Yansong Feng, Yu Hong, and Ruifang He, Eds., Cham, 2021, pp. 155–168, Springer International Publishing.
  * [21] Zhexi Zhang, Wei Zhu, Junchi Yan, Peng Gao, and Guowang Xie,  “Automatic student network search for knowledge distillation,”  2020 25th International Conference on Pattern Recognition (ICPR), pp. 2446–2453, 2021.
  * [22] Wei Zhu, Yuan Ni, Xiaoling Wang, and Guo Tong Xie,  “Discovering better model architectures for medical query understanding,”  in North American Chapter of the Association for Computational Linguistics, 2021.
  * [23] Mostafa Dehghani et al.,  “Universal transformers,”  arXiv, vol. abs/1807.03819, 2018.
  * [24] Zhen Zhang, Wei Zhu, Jinfan Zhang, Peng Wang, Rize Jin, and Tae-Sun Chung,  “Pcee-bert: Accelerating bert inference via patient and confident early exiting,”  in NAACL-HLT, 2022.
  * [25] Teerapittayanon et al.,  “Branchynet: Fast inference via early exiting from deep neural networks,”  in 2016 23rd International Conference on Pattern Recognition (ICPR), 2016, pp. 2464–2469.
  * [26] Weijie Liu et al.,  “Fastbert: a self-distilling BERT with adaptive inference time,”  arXiv, vol. abs/2004.02178, 2020.
  * [27] Ji Xin et al.,  “Deebert: Dynamic early exiting for accelerating BERT inference,”  arXiv, vol. abs/2004.12993, 2020.
  * [28] Yigitcan Kaya et al.,  “How to stop off-the-shelf deep neural networks from overthinking,”  CoRR, vol. abs/1810.07052, 2018.
  * [29] Roy Schwartz et al.,  “The right tool for the job: Matching model and instance complexities,”  arXiv, vol. abs/2004.07453, 2020.
  * [30] Ji et al. Xin,  “BERxiT: Early exiting for BERT with better fine-tuning and extension to regression,”  in Proceedings of the 16th Conference of the European Chapter of the Association for Computational Linguistics: Main Volume, Online, Apr. 2021, pp. 91–104, Association for Computational Linguistics.
  * [31] Tal Schuster et al.,  “Consistent accelerated inference via confident adaptive transformers,”  CoRR, vol. abs/2104.08803, 2021.
  * [32] Wei Zhu,  “LeeBERT: Learned early exit for BERT with cross-level optimization,”  in Proceedings of the 59th Annual Meeting of the Association for Computational Linguistics and the 11th International Joint Conference on Natural Language Processing (Volume 1: Long Papers), Online, Aug. 2021, pp. 2968–2980, Association for Computational Linguistics.
  * [33] Wei Zhu, Xiaoling Wang, Yuan Ni, and Guotong Xie,  “GAML-BERT: Improving BERT early exiting by gradient aligned mutual learning,”  in Proceedings of the 2021 Conference on Empirical Methods in Natural Language Processing, Online and Punta Cana, Dominican Republic, Nov. 2021, pp. 3033–3044, Association for Computational Linguistics.
  * [34] Geoffrey E. Hinton et al.,  “Distilling the knowledge in a neural network,”  ArXiv, vol. abs/1503.02531, 2015.
  * [35] Alex Wang et al.,  “GLUE: A multi-task benchmark and analysis platform for natural language understanding,”  CoRR, vol. abs/1804.07461, 2018.
  * [36] Alice Coucke et al.,  “Snips voice platform: an embedded spoken language understanding system for private-by-design voice interfaces,”  arXiv, vol. abs/1805.10190, 2018.
  * [37] Hemphill et al.,  “The ATIS spoken language systems pilot corpus,”  in Speech and Natural Language: Proceedings of a Workshop Held at Hidden Valley, Pennsylvania, June 24-27,1990, 1990.
  * [38] Pengcheng Yang et al.,  “SGM: sequence generation model for multi-label classification,”  arXiv, vol. abs/1806.04822, 2018.
  * [39] Jeff Atwood,  “Stack overflow creative commons data dump,”  2009\.
  * [40] Haoli Bai, Wei Zhang, Lu Hou, Lifeng Shang, Jing Jin, Xin Jiang, Qun Liu, Michael R. Lyu, and Irwin King,  “Binarybert: Pushing the limit of bert quantization,”  in Annual Meeting of the Association for Computational Linguistics, 2020.
  * [41] Ilya Loshchilov et al.,  “Decoupled weight decay regularization,”  in ICLR, 2019.
  * [42] Wolf et al.,  “Transformers: State-of-the-art natural language processing,”  in Proceedings of the 2020 conference on empirical methods in natural language processing: system demonstrations, 2020, pp. 38–45.