# Single–pulse studies of three millisecond pulsars

N. T. Palliyaguru,1 B. B. P. Perera,2 M. A. McLaughlin,3,4 S. Osłowski,5 and
G. L. Siebert6

1 Department of Physics and Astronomy, Texas Tech University, Lubbock, TX
79409, USA
2 Arecibo Observatory, University of Central Florida, HC3 Box 53995, Arecibo,
PR 00612, USA
3 Department of Physics and Astronomy, West Virginia University, Morgantown,
WV 26501, USA
4 Center for Gravitational Waves and Cosmology, West Virginia University,
Chestnut Ridge Research Building, Morgantown, WV 26505, USA
5 Manly Astrophysics, 15/41-42 East Esplanade, Manly, NSW 2095, Australia
6 Department of Physics, University of Wisconsin Madison, Madison, WI 53703,
USA E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

Single–pulse studies are important to understand the pulsar emission mechanism
and the noise floor in precision timing. We study total intensity and
polarimetry properties of three bright millisecond pulsars – PSRs J1022+1001,
J1713+0747, and B1855+09 – that have detectable single pulses at multiple
frequencies. We report for the first time the detection of single pulses from
PSRs J1022+1001 and J1713+0747 at 4.5 GHz. In addition, for those two pulsars
the fraction of linear polarization in the average profile is significantly
reduced at 4.5 GHz, compared to 1.38 GHz, which could support the expected
deviation from a dipolar field closer to the pulsar surface. There is a hint
of orthogonal modes in the single pulses of PSR J1713+0747. More sensitive
multi–frequency observations may be useful to confirm these findings. The
jitter noise contributions at 1.38 GHz, scaled to one hour, for PSRs
J1022+1001, J1713+0747 and PSR B1855+09 are $\approx$135 ns, $\approx$45 ns,
and $\approx$60 ns respectively and are consistent with previous studies. We
also show that selective bright-pulse timing of PSR J1022+1001 yields improved
root-mean-square residuals of $\approx$22 $\mu$s, which is a factor of
$\approx$3 better than timing using single pulses alone.

###### keywords:

stars: neutron – pulsars: general – pulsars: individual: PSR J1022+1001,
J1713+0747, B1855+09

††pubyear: 2021††pagerange: Single–pulse studies of three millisecond
pulsars–LABEL:lastpage

##  1 Introduction

In general, pulsar emission remains poorly understood. The commonly accepted
model suggests that emission is produced by curvature radiation from bunches
of charged particles moving along magnetic field lines at relativistic speeds.
The emitted photons split into electron–positron pairs, producing a cascade of
secondary pair plasma (Goldreich & Julian, 1969; Ruderman & Sutherland, 1975).
Components of the electric field that are parallel ($\rm E_{\parallel}$) and
perpendicular ($\rm E_{\perp}$) to the magnetic field lines produce highly
linearly polarized emission. Circular polarization is thought to be generated
both by intrinsic mechanisms and/or by propagation effects within the
magnetosphere (Melrose, 1995), where a time delay is introduced between $\rm
E_{\parallel}$ and $\rm E_{\perp}$ due to different refractive indices of the
two orthogonal components in the magnetosphere (Gil & Snakowski, 1990). As the
pulsar beam crosses our line of sight, sweeping through different magnetic
field lines, the polarization position angle (PPA), which is the angle between
the magnetic field line and the fiducial plane (i.e., the plane passing
through the rotation and magnetic axes) makes an S–shaped curve (Radhakrishnan
& Cooke, 1969). The PPA information is utilized in constraining the geometry
of pulsars using the rotating vector model (e.g. Rankin, 1983a, b; Mitra &
Rankin, 2011).

Millisecond pulsars (MSPs) are old neutron stars that are spun up to
millisecond periods. They have smaller magnetospheres and relatively smaller
($\sim$four orders of magnitude) magnetic field strengths than canonical
pulsars, which have not been through the recycling process (Lorimer & Kramer,
2012). This poses the question of whether emission of MSPs is different from
canonical pulsars. For example, the emission beam size of a canonical pulsar
is correlated with its period but such a correlation cannot be clearly seen
for MSPs (Kramer et al., 1998). In addition, MSP profiles appear to show less
evolution with frequency than those of canonical pulsars (Kramer et al.,
1998). Some properties of canonical and millisecond pulsars are similar,
however. For instance, MSPs can emit giant pulses (Heiles & Campbell, 1970),
like normal pulsars; there is evidence that giant pulses are linked to high
magnetic-field strengths at the pulsar light cylinder and to high-energy
emission (Romani & Johnston, 2001; Enoto et al., 2021). MSPs too may show a
deviation from the dipolar field line structure close to the NS surface (e.g.
Kramer et al., 1998; Gil et al., 2002; Kalapotharakos et al., 2021).

Pulsar emission is extremely complex, requiring sophisticated models to
explain observations that do not fit into the general picture. Average pulse
profiles are stable in general, with some notable exceptions, such as PSRs
J1713+0747 (Brook et al., 2018; Xu et al., 2021), B1937+21 (Brook et al.,
2018), and PSR B1828-11 (Lyne et al., 2010). However, emission is highly
variable at a single–pulse level. Phenomena such as giant pulses (Heiles &
Campbell, 1970), mode changing (Backer, 1970), drifting sub-pulses (Drake &
Craft, 1968; Backer, 1973), nulling (Rankin, 1986) and microstructure within
single pulses (Cordes et al., 1990) are not easily explained by a simple
emission model. The PPAs of many pulsars, in particular MSPs, deviate from the
S–shaped sweep and show discontinuities or unexpected jumps. These are
sometimes separated by $90^{\circ}$, named orthogonal polarization modes
(OPMs), and may result from either two highly polarized orthogonal modes of
emission (e.g. Backer et al., 1976; Cordes et al., 1978) or propagation
effects (McKinnon & Stinebring, 2000). Even though the curvature radiation
with bunches explains the radio emission of pulsars to some extent, it has
serious drawbacks, including the inability to explain the production of
bunches and their existence over the characteristic time of emission (Melrose,
1992). Furthermore, secondary pair plasma is insufficient to reach the charge
density required by Maxwell’s equations (e.g. Melrose, 1995). As alternative
mechanisms, relativistic plasma emission (Melrose & Gedalin, 1999), anomalous
Doppler emission (e.g. Machabeli & Usov, 1979), linear acceleration emission
(Cocke, 1973), and maser emission (Zhelezniakov & Shaposhnikov, 1979) have
been proposed (see Melrose et al., 2021, for a review of these effects).
However, none of these theoretical methods are capable of fully explaining the
pulsar emission mechanism (Melrose, 1995; Melrose et al., 2021). Single-pulse
studies can provide valuable diagnostics of emission mechanisms. For example,
Crab pulsar microstructures suggest coherent emission from strong plasma
turbulence (Hankins & Eilek, 2007). Also, highly polarized single pulses have
been used to distinguish between coherent curvature radiation and maser
emission in the past (Mitra et al., 2009). Therefore, single–pulse studies,
including single–pulse–polarimetry, of pulsars with a variety of properties
are crucial to fully understand pulsar emission physics (Rankin, 1986).

MSPs are also important for low-frequency, from nHz to $\mu$Hz, gravitational
wave detection experiments through pulsar timing arrays (PTAs), which require
precise measurements of pulse times of arrival (TOAs) and sub-microsecond
timing accuracy (Alam et al., 2021; Desvignes et al., 2016; Kerr et al., 2020;
Perera et al., 2019, 2018). Pulse–to–pulse jitter, profile variations,
variation in polarization properties, and polarization calibration errors are
among phenomena that could affect the precision of TOAs (Cordes & Shannon,
2010; Liu et al., 2011; Osłowski et al., 2013). These can contribute to noise
that can limit the timing precision of MSPs and therefore sensitivity of PTAs
to gravitational waves. Therefore studying MSPs at a single–pulse level to
understand the contributions of these various phenomena can offer useful
insights.

Single–pulse studies of MSPs have been sparse due to signal–to–noise ratio
limitations due to their low fluxes (Kramer et al., 1998) and data acquisition
requirements such as the need for high time resolution sampling. The few
previous studies have revealed highly linearly polarized single pulses and
sub–pulse microstructure in J0437$-$4715 (Jenet et al., 1998; Osłowski et al.,
2014; De et al., 2016), giant pulses from B1937$+$21 (McKee et al., 2019),
pulse jittering from PSR J1713$+$0747 (Shannon & Cordes, 2012) and PSR
J1022+1001 (Liu et al., 2015; Feng et al., 2020), and sub-pulse drifting from
PSR J1713$+$0747 (Liu et al., 2015).

Polarization information in single pulses could get lost when averaging,
causing depolarization in averaged pulse profiles (Backer, & Rankin, 1980).
Therefore studying single pulses may provide clues to the emission physics of
pulsars. In addition to emission mechanisms, there are other advantages to
performing polarimetry. Britton (2000) suggested the possibility of using
invariant profiles to avoid errors due to calibration and van Straten (2006)
suggested the use of polarimetric profiles to improve timing. Furthermore,
Osłowski et al. (2013) used polarization information to correct for
pulse–to–pulse variability in PSR J0437$-$4715, which resulted in a 40%
improvement in the timing precision of the pulsar. In light of these ideas, we
investigate single–pulse emission properties in three stable millisecond
pulsars – PSRs J1022+1001, J1713+0747, and B1855+09 – which are monitored
regularly by PTAs (see Perera et al., 2019). PSR J1022+1001 is known to show
long–term profile instabilities (e.g. Kramer et al., 1999; Liu et al., 2015;
Padmanabh et al., 2021) and PSR J1713+0747 recently underwent a significant
pulse shape change (e.g. Xu et al., 2021). While the origin of such changes is
likely intrinsic, and not due to incorrect calibration, these events can
introduce variations in timing residuals, degrading the PTA sensitivity to
gravitational waves.

The paper is structured as follows. The details of observations and data
processing are presented in Section 2. Data analysis techniques are presented
in Section 3 and results are presented in Section 4. The conclusion and
discussion of the study are presented in Section 5.

##  2 Observations and data processing

PSRs J1022$+$1001, J1713+0747, and B1855+09 were observed with the 305-m
William E. Gordon Telescope at the Arecibo observatory in Puerto Rico. The
basic parameters of these three millisecond
pulsars111https://www.atnf.csiro.au/research/pulsar/psrcat/ are given in Table
LABEL:tb:msplist. Arecibo observations were carried out between August 5, 2018
and August 25, 2019. Data were recorded using the Puerto Rico Ultimate Pulsar
Processing Instrument (PUPPI) at center frequencies of 430 MHz, 1380 MHz, 2030
MHz and 4500 MHz in full Stokes mode with 8–bit sampling. Our multi-frequency
observation details, including the usable bandwidth, number of channels across
the usable bandwidth and the sampling time, for each pulsar, are listed in
Table LABEL:tb:obs. The usable bandwidth at 1.38 GHz and 2.3 GHz is less than
the full receiver bandwidth due to radio frequency interference (RFI). While
baseband recording mode allows higher time resolution, we accumulated data in
search mode in order to maximize the available bandwidth and thereby to
improve the S/N of the pulsar data. PSRs J1713+0747 and B1855+09 data were
coherently dedispersed using dspsr222http://dspsr.sourceforge.net (van Straten
& Bailes, 2011) at the dispersion measure listed in Table LABEL:tb:obs. PSR
J1022+1001 was observed in the incoherent search mode. Each observation
session started with a noise calibrator injection scan followed by the pulsar
observation.

Table 1: Properties of the observed pulsars: The pulsar period, period
derivative, the dispersion measure, and references.

Name | Period | $\dot{P}$ | DM | Reference
---|---|---|---|---
| (ms) | ($\rm s\,s^{-1}$) |  $(\rm cm^{-3}\,pc$) |
J1022+1001 | 16.453 | $4.334\times 10^{-20}$ | 10.252 | 1,2
J1713+0747 | 4.570 | $8.530\times 10^{-21}$ | 15.917 | 3,4
B1855+09 | 5.362 | $1.784\times 10^{-20}$ | 13.314 | 4

  * •

References: (1) Hotan et al. (2006), (2) Reardon et al. (2016), (3) Zhu et al.
(2015), (4) Arzoumanian et al. (2018).

The recorded data were processed with dspsr (van Straten & Bailes, 2011) to
obtain single pulses, which were further processed with
psrchive333http://psrchive.sourceforge.net routines (Hotan et al., 2004). The
data from each pulsar were processed with 512 pulse phase bins given the time
resolution used in the observation setup. The automatic median zapping
algorithms of psrchive were used to remove narrow–band and impulsive RFI.
Remaining RFI was removed by visual inspection.

Polarization calibration was performed to correct for the differential gain
and phase between the two polarization channels of the receiver due to
imperfections in the amplifiers and mismatch between the cable chains along
the two paths, assuming an ideal feed. For this, we used the calibrator scan
which is a 25 Hz winking cal signal injecting a fully linearly polarized
signal in between the two polarization probes at a 45∘ position angle. Next,
the instrumental response was determined using polarization calibration
modeling (pcm) as outlined in Britton (2000), assuming equal ellipticities of
receptors. A flux calibrator was not observed during these observations. The
polarization properties of the integrated profiles are shown in Figure 1,
which are consistent with previously published profiles (Dai et al., 2015).

Faraday rotation due to the magnetized ISM and the ionosphere causes a change
in the PPA with frequency (e.g. Simard-Normandin & Kronberg, 1980; Lyne &
Smith, 1989). The change in the PPA due to Faraday rotation is given by

$\Delta\Psi=\lambda^{2}\times\rm RM,$ (1)

where $\lambda$ is the observation wavelength. The rotation measure (RM),
which depends on the average magnetic field and the electron density along the
line–of–sight, was calculated using the rmfit program in psrchive, which
searches for a peak in the linear polarization for trial RMs. The best-fit RM
values (published values in parenthesis) are given in Table LABEL:tb:obs.
Given the $\lambda^{2}$ dependence in Equation 1, and the decrease in flux
density with increasing frequency, the sensitivity to the change in the PPA
decreases with increasing frequency. Therefore, the RM values derived from
1.38 GHz data were used to correct the data at higher frequencies.

Our RM values at 1.38 GHz are consistent within errors with published values
for MSPs J1713+0747 and B1855+09 (see Table LABEL:tb:obs). However,
significant trends in the measured RM has been observed for PSR J1713+0747
previously (Wahl et al., 2021). Our measured RM for PSR J1022+1001 differs
significantly from its published value, however, PSR J1022+1001 is known to
show changes in the RM when passing close to the sun i.e. when the angular
separation between the pulsar and the sun is $<3^{\circ}$ (You et al., 2012)
and also shows significant long–term RM variations (Yan et al., 2011). The
angular separation between the sun and the pulsar was $\approx 15^{\circ}$
during the 1.38 GHz observations. While it is not possible to pinpoint to the
reason for the discrepancy between our measured and previously published RM
values for PSR J1022+1001, we note that Noutsos et al. (2015), Feng et al.
(2020) and Dai et al. (2015) also find higher RM values of 2.18(2), 2.9(2) and
4.68(6) $\rm rad\,m^{-2}$ respectively.

For the pulse phase jitter analysis, average profiles of 50, 100, 200, 500,
and 1000 pulses were created. A TOA for each profile was generated by cross-
correlating it against a noise-free template profile. This template was
obtained by fitting Gaussian components to a high S/N profile obtained by
averaging over the full observation for a given frequency. Timing residuals
were obtained using the tempo2 pulsar timing package (Edwards et al., 2006;
Hobbs et al., 2006). The timing ephemerides are obtained from data published
in Perera et al. (2019).

Table 2: Observation information, including MJD, frequency, observation
length, bandwidth, number of channels across the bandwidth, sampling time,
receiver temperature, sky temperature, measured rotation measures (published
rotation measures are given in parentheses), equivalent width, measured mean
flux density from the radiometer equation, and the number (and fraction) of
pulses detected ($\rm S/N>5$).

Name | MJD | $\nu$ | Length | BW | $\rm N_{chan}$ | G | $t_{\rm samp}$ |  Trec |  Tsky | RM |  Weq | $S_{\rm mean}$ | N
---|---|---|---|---|---|---|---|---|---|---|---|---|---
|  | (MHz) | (s) | (MHz) |  | K/Jy |  $(\mu$s) | (K) | (K) | $(\rm rad\,m^{-2})$ | (ms) | (mJy) |
J1022+1001 | 58335 | 1380 | 300 | 600 | 384 | 8 | 20.48 | 30 | 0.62 |  $5.31\pm 0.04$ (–0.3, 2.18) |  1.028$\pm$0.002 |  4.89$\pm$0.15 |  10755 ($\approx 59\%$)
| 58651 | 4500 | 1800 | 800 | 512 | 4 | 40.96 | 30 | 0.03 | – |  0.577$\pm$0.003 |  0.90$\pm$0.11 |  1175 ($\approx 1\%$)
J1713+0747 | 58697 | 1380 | 300 | 600 | 384 | 8 | 10.24 | 30 | 2.92 |  $10.69\pm 0.04$ (13$\pm$ 2) |  0.1752$\pm$0.0003 |  10.18$\pm$0.17 |  64458 ($\approx 98\%$)
| 58697 | 4500 | 300 | 800 | 512 | 4 | 10.24 | 30 | 0.14 | – |  0.195$\pm$0.002 |  1.37$\pm$0.27 |  53 ($\approx 0.08\%$)
B1855+09 | 58720 | 430 | 600 | 20 | 64 | 11 | 2.56 | 50 | 109.9 | – |  0.53$\pm$0.01 |  10.61$\pm$2.15 | –
| 58720 | 1380 | 1200 | 600 | 384 | 8 | 10.24 | 30 | 5.11 |  $24.99\pm 0.18$ (20$\pm$4) |  0.483$\pm$0.001 |  5.21$\pm$0.09 |  11480 ($\approx 10\%$)
| 58720 | 2030 | 600 | 460 | 384 | 8 | 10.24 | 40 | 1.94 | – |  0.43$\pm$0.02 |  0.39$\pm$0.17 | –

  * •

Notes:

  * •

The published RM values are obtained from Yan et al. (2011), Noutsos et al.
(2015) and Gentile et al. (2018).

  * •

The RM of PSR B1855+09 is not calculated at 430 MHz due to the small
bandwidth.

  * •

$\rm T_{\rm rec}$ and Gain values for AO are obtained from
http://www.naic.edu/~astro/RXstatus/rcvrtabz.shtml.

  * •

$\rm T_{\rm sky}$ in the direction of the pulsar is calculated according to
Haslam et al. (1982).

##  3 Data analysis

In this section we discuss the data analysis tools used to measure total
intensity properties, polarimetry of single pulses, and pulse phase jitter of
PSRs J1022+1001, J1713+0747, and B1855+09.

Table LABEL:tb:obs lists the measured equivalent pulse width $W_{eq}$ and the
measured flux density $S_{mean}$ of the average profile over the full
observation length from our multi-frequency observations. $W_{\rm eq}$ is
defined as the width of a top hat pulse with the peak amplitude $I_{\rm peak}$
and the same area as the on–pulse region and is calculated as

$W_{\rm eq}=\frac{\sum_{i=n_{1}}^{N}I_{i}}{I_{\rm peak}},$ (2)

where $I_{i}$ is the intensity of the $i^{th}$ bin of the on–pulse region
ranging from bins $n_{1}$ to N in baseline-corrected data. The error of
$W_{\rm eq}$ is calculated by applying error propagation on Equation 2 and
using the off and on–pulse rms as the errors of $I_{i}$ and $I_{\rm peak}$
respectively.

From the radiometer equation, the root mean square noise fluctuation is given
by

$\Delta S_{\rm sys}=\frac{T_{\rm
sys}}{\,G\sqrt{N_{p}\,t_{obs}\,\Delta\nu}}=\rm C\,\sigma_{p}.$ (3)

Here, C is the scaling factor, $\sigma_{p}$ is the standard deviation of the
off–pulse region, $T_{\rm sys}$ is the system temperature, which is the sum of
the sky temperature ($T_{\rm sky}$), the receiver temperature ($T_{\rm rec}$),
and the telescope spill over, $G$ is the telescope gain, $\Delta\nu$ is the
bandwidth, $t_{\rm obs}$ is the observation length, and $N_{p}=2$ is the
number of polarization channels (see Lorimer & Kramer, 2012). The measured
mean flux density $S_{\rm mean}$ is calculated by scaling the profile by $\rm
C$ (i.e. $\Delta S_{\rm sys}/\sigma_{p})$. The values for $T_{\rm rec}$,
$T_{\rm sky}$, and $G$ for different frequencies are listed in Table
LABEL:tb:obs.

###  3.1 Single–pulse properties

Figure 1:  Multi-frequency single pulse observations of PSRs J1022$+$1001,
J1713$+$0747, and B1855$+$09 at observing frequency 1.38 and 4.5 GHz (see
different columns). For each pulsar at each frequency, a sequence of 2000
single pulses in grey scale (top panel) and the averaged profiles of each
pulse sequence, including total intensity (black), linear (red) and circular
(blue) polarization (second panel from top) are shown. Three highest S/N
single-pulse profiles are shown from each observation (bottom three panels). A
single pulse associated with the leading component of PSR J1022+1001 at 1.38
GHz with a low fraction of linear polarization is also shown. Two single
pulses of PSR J1022+1001 at 4.5 GHz which appear at the small trailing
component at phase $\approx 0.53$ (away from the main component) are also
shown.

We analyze single–pulse properties of the five data sets. Single-pulse
amplitude and S/N distributions are used to identify phenomena such as giant
pulses, which have flux density $>10\times$ the mean flux density (Karuppusamy
et al., 2012; Knight, 2007). The S/N for each single-pulse is calculated as

$S/N=\frac{I_{\rm peak}}{\sigma_{p}}$ (4)

where $I_{\rm peak}$ is the peak amplitude, and $\sigma_{p}$ is the off–pulse
standard deviation of the pulse profile. The number and the fraction of
detected single pulses ($\rm S/N>5$) is listed in Table LABEL:tb:obs.

Figure 2:  Histograms of peak amplitude (top) and S/N (bottom) of PSR
J1022+1001 at 1.38 GHz (a) and 4.5 GHz (b), PSR J1713+0747 at 1.38 GHz (c),
and 4.5 GHz (d), and PSR B1855+09 at 1.38 GHz (e). For PSR J1022+1001 at 1.38
GHz, the histograms of single–pulses associated with the leading (black) and
trailing (blue) components are also shown separately in addition to the
histograms of all single pulses (red). The brightest seem to correspond to the
trailing component.

Figure 3:  Measured single–pulse equivalent width and S/N values for PSR
J1022+1001 at 1.38 GHz (a) and 4.5 GHz (b), PSR J1713+0747 at 1.38 GHz (c),
and PSR B1855+09 at 1.38 GHz (d). The S/N cutoff is used to show only
measurements > $3\sigma$. Widths of PSR J1713+0747 single pulses at 4.5 GHz
are not shown due to the small fraction of detected pulses with S/N$>$5\.

###  3.2 Single–pulse polarimetry

Stokes parameters I, Q, U, V describe the polarization of electromagnetic
waves, such that I is the total intensity, Q and U describe the intensity of
linear polarization L as $L=\sqrt{Q^{2}+U^{2}}$, and V is the intensity of
circular polarization (see Lorimer & Kramer, 2012). The polarization position
angle is defined as $\rm PPA=0.5\tan^{-1}(U/Q)$.

Table 3: Polarization parameters: percentage fractional linear and circular
polarization of the average profile and the peak of the histograms of linear
and circular polarizations shown in Figure 4.

Pulsar | Frequency | $\rm L_{\rm av}/I_{\rm av}$ | $\rm V_{\rm av}/I_{\rm av}$ | $\rm L/I$ | $\rm V/I$
---|---|---|---|---|---
| (MHz) |  |  |  |
J1022+1001 | 1380 |  73.9$\pm$0.2 |  –18.0$\pm$0.2 |  25.9$\pm$1.5 , 73.8$\pm$1.5 |  –14.6$\pm$1.6
J1022+1001 | 4500 |  3.2$\pm$0.3 |  6.9$\pm$0.4 |  26.5$\pm$0.7 |  6.0$\pm$1.3
J1713+0747 | 1380 |  36.6$\pm$0.5 |  0.9$\pm$0.5 |  36.3$\pm$0.9 |  –0.5$\pm$1.5
J1713+0747 | 4500 |  30.1$\pm$0.6 |  –1.4$\pm$0.8 |  34.5$\pm$0.8 |  –0.6$\pm$1.2
B1855+09 | 1380 |  6.6$\pm$0.1 |  –3.5$\pm$0.2 |  27.3$\pm$0.9 |  –4.6$\pm$1.7

  * •

For PSR J1022+1001 at 1.38 GHz the peak fractional linear polarization for
leading and trailing components are listed separately.

Figure 4:  Histograms of linear polarization (top) and circular polarization
(bottom) for PSR J1022+1001 at 1.38 GHz (a) and 4.5 GHz (b), PSR J1713+0747 at
1.38 GHz (c), and 4.5 GHz (d), and PSR B1855+09 at 1.38 GHz (e). For PSR
J1022+1001 at 1.38 GHz, the histograms of single–pulses associated with the
leading (black) and trailing (blue) components are also shown separately.
Single pulses associated with the trailing component show a higher fractional
linear polarization than those associated with the leading component.

Figure 5:  The total intensity average profile of the full dataset (top) and
the single-pulse PPA distribution (bottom grey scale) of for PSR J1022+1001 at
1.38 GHz (a) and 4.5 GHz (b), PSR J1713+0747 at 1.38 GHz (c), and 4.5 GHz (d),
and PSR B1855+09 at 1.38 GHz (e). The PPA of the average profile (orange) is
also marked.

###  3.3 Pulse phase jitter

Figure 6:  Correlation between TOAs of top and bottom halves of the band for
1000-pulse averages along with the best-fit line (top), the cross correlation
function (middle), and the rms timing residuals from jitter vs the number of
pulses averaged (bottom) for PSR J1022+1001 at 1.38 GHz (a) and 4.5 GHz (b),
PSR J1713+0747 at 1.38 GHz (c), and 4.5 GHz (d), and PSR B1855+09 at 1.38 GHz
(e). In the bottom plots the best-fit line to jitter noise is shown by the
dashed line. The correlation is negligible for both PSRs J1022+1001 and
J1713+0747 at 4.5 GHz for 1000-pulse averages likely due to radiometer noise
in profiles.  Figure 7:  The rms timing residuals of PSR J1022+1001
single–pulses of $\rm S/N>25$ at 1.38 GHz.

As demonstrated by Shannon & Cordes (2012), timing residuals due to
pulse–to–pulse jitter are expected to be correlated across radio frequency.
For average profiles formed with $N$ pulses, the root-mean-square (rms)
residuals due to jitter noise, $\sigma_{J}(N)$, is quantified as the
quadrature difference between the measured rms timing residuals, $\sigma_{\rm
obs}(N)$, and the expected rms timing residuals from radiometer noise only,
$\sigma_{\rm rad}(N)$, such that (Shannon & Cordes, 2012)

$\sigma_{J}(N)^{2}=\sigma_{\rm obs}(N)^{2}-\sigma_{\rm rad}(N)^{2}.$ (5)

In order to determine the effect of jitter, composite profiles were formed
with $N=$ 50, 100, 200, 500, and 1000 pulses averaged together. For these
various $N$ values, $\sigma_{\rm obs}(N)$ was calculated from the timing
residuals of the measured TOAs (obtained as described in Section 2) and
$\sigma_{\rm rad}(N)$ was calculated from simulated profiles of similar S/N
levels as the composite profiles (Downs & Reichley, 1983). The simulated
profiles were constructed by adding white noise to a noise-free template with
signal-to-noise drawn from the signal-to-noise distribution of the observed
composite profiles. We note that the simulated dataset does not contain
uncertainties due to profile shape variations or any other phenomena. The rms
residuals due to jitter noise were determined from Equation 5 for various $N$
values.

Timing residuals due to jitter noise are correlated between sub-bands. The
band may be split into several sub-bands and the correlation may be calculated
between adjacent sub-bands (Parthasarathy et al., 2021). Following Shannon &
Cordes (2012), we calculate the rms due to jitter by correlating timing
residuals from the top and bottom halves of the band as
$\sigma_{j}=\sqrt{\rm{CCF(0)}}$. Here, CCF(0) is the zero–lag value of the
cross–correlation function. A different template profile was used for each
subband.

Table 4: Phase jitter analysis: Pulsar name, frequency, the slope of the best fit to the timing residual correlation for 1000–pulse averages ($m$ -see the top panel of Figure 6), the corresponding zero–lag value of the CCF (CCF(0)), rms residuals due to jitter from the CCF ($\sigma_{\rm j,CCF}$), measured rms timing residuals ($\sigma_{\rm obs}$), the reduced chi-square of the timing residuals ($\chi_{r}^{2}$), rms residuals due to radiometer noise ($\sigma_{\rm rad}$), rms residuals due to jitter ($\sigma_{\rm j}$) from the quadrature difference between $\sigma_{\rm obs}$ and $\sigma_{\rm rad}$, jitter scaling parameters (A, $\beta$ – see text), and $\sigma_{\rm j}$ scaled to one hour ($\sigma_{\rm j,1\,hr}$).  Pulsar | Frequency | $m$ | CCF(0) | $\sigma_{\rm j,CCF}$ | $\sigma_{\rm obs}$ | $\chi_{r}^{2}$ | $\sigma_{\rm rad}$ | $\sigma_{\rm j}$ |  (A, $\beta$) | $\sigma_{\rm j,1\,hr}$
---|---|---|---|---|---|---|---|---|---|---
| (MHz) |  |  ($\mu$s2) | (ns) | (ns) |  | (ns) | (ns) |  | (ns)
J1022+1001 | 1380 | 0.71 | 4.7 | 2164 | 2086 | 16.1 | 467 | 2033 |  (61.10,$-$0.50) |  135$\pm 29$
J1022+1001 | 4500 | 0.18 | 25.3 | 5033 | 4502 | 1.69 | 3764 | 2469 |  (146.4,$-0.58$) |  117$\pm$52
J1713+0747 | 1380 | 1.04 | 0.72 | 846 | 851 | 49.9 | 121 | 842 |  (17.03,$-$0.44) |  45$\pm$7
J1713+0747 | 4500 | 0.03 | 0.22 | 473 | 1939 | 1.30 | 1779 | 772 |  (43.15,$-$0.57) |  17$\pm$22
B1855+09 | 1380 | 0.80 | 3.31 | 1819 | 2003 | 7.87 | 914 | 1782 |  (70.39,$-$0.53) |  60$\pm$25

##  4 Results

In this section we first present overall single–pulse total intensity,
polarimetry, and jitter results and then present the results for individual
pulsars in subsections.

Single–pulse properties: Single pulses of $\rm S/N>5$ are detected for PSRs
J1022+1001 and J1713+0747 at 1.38 GHz and 4.5 GHz, and for B1855+09 at 1.38
GHz. The 2.3 GHz data set of B1855+09 does not show detectable single pulses.
To the best of our knowledge, this is the first report of the detection of
single pulses at frequencies $>2$ GHz in MSPs.

Figure 1 shows a single-pulse sequence (a stack of individual pulses in time
vs pulse–phase), the total intensity and polarization profiles (see Section
3.2 for polarimetry) of averaged data and the three highest S/N individual
single pulses from each data set. Figure 2 shows the histograms of peak pulse
amplitude distribution and the S/N distribution of single pulses of PSRs
J1022+1001, J1713$+$0747, and B1855+09 at frequencies 1.38 GHz and 4.5 GHz.

Figure 3 shows the single–pulse equivalent width vs S/N. Overall, single
pulses with higher S/N seem to be narrower than the ones with smaller S/N at
1.38 GHz, as shown in Figure 3. In general, Gaussians or a series of Gaussians
are fit to pulse profiles to estimate the width. The equivalent width method
is better for automated calculation for many pulses, where goodness of fit is
difficult to ascertain for each pulse, and to facilitate width measurements
for weak pulses. At low S/N, the true width may be difficult to measure
resulting in large width estimates. In Figure 3, we plot width measurements
determined with > $3\sigma$ significance.

Polarimetry: Figure 4 shows the histograms of single–pulse linear and circular
polarization fraction at the peak. The bias in L has been removed according to
Everett & Weisberg (2001). Table LABEL:tb:pol lists the linear and circular
polarization fraction of the average profiles and the peak of the
distributions of single-pulse polarization fractions shown in Figure 4. The
error on the polarization fraction of the average profile are obtained by
using the standard deviation of I, L, and V across the profile and applying
error propagation.

Phase resolved histograms of the PPA for single–pulse data were obtained using
psrspa routine in psrchive. Figure 5 shows the logarithm distribution of
single-pulse PPA as a function of pulse phase along with the PPA of the
average profile. The histogram is weighted by the total intensity.

Pulse phase jitter: Figure 6 shows the correlation between the timing
residuals from the top and bottom halves of the band, for 1000-pulse averages,
along with the best-fit line whose slope should be ideally be unity. The
slopes of the best-fit lines for the three pulsars are listed in Table
LABEL:tb:jitter. Table LABEL:tb:jitter also lists the zero lag value of the
cross–correlation function CCF(0), rms timing residuals $\sigma_{\rm obs}$ for
1000-pulse averages, the reduced chi-square of the timing residuals
($\chi_{r}^{2}$), the TOA error from radiometer noise $\sigma_{\rm rad}$, rms
residuals due to jitter calculated from Equation 5, the best–fit parameters of
the line $\sigma_{J}(N)$ vs $N$, and rms residuals due to jitter scaled to one
hour. The rms residuals due to jitter for one hour averages were determined
from the parameters of the best–fit line to $\sigma_{J}(N)$ vs $N$ such that
$\sigma_{J}(N)=A\,N^{\beta}$, as shown in the bottom panel of Figure 6. The
error on the jitter level scaled to one hour is obtained by applying error
propagation to $\sigma_{J}(N)=A\,N^{\beta}$. The quoted error bars correspond
to $1\sigma$ level.

In the quadrature difference method, measuring the radiometer noise component
accurately is important to get a reliable estimate for jitter noise since the
datasets at higher frequencies are in general of low S/N. Higher S/N for these
pulsars would not be achievable with any other existing telescope at 4.5 GHz.
A comparison of the jitter level from the cross-correlation method and
quadrature difference ($\sigma_{\rm j,CCF}$ and $\sigma_{j}$ in Table
LABEL:tb:jitter) between the observed and simulated rms timing residuals shows
that the values are comparable (within a factor of $\lesssim 1.1$) in all
datasets at 1.38 GHz. The inconsistencies at 4.5 GHz in both datasets of
J1022+1001 and J1713+0747 could perhaps be because the radiometer noise
component is underestimated or because jitter is decorrelated across the band
(e.g. Parthasarathy et al., 2021).

###  4.1 PSR J1022+1001

As shown in Figure 1, PSR J1022+1001 shows a double–peaked profile with a
highly linearly polarized ($\approx 74\%$) trailing component at 1.38 GHz, as
also noted previously (Hotan et al., 2004; Liu et al., 2015). At 4.5 GHz the
leading component has reduced in amplitude with a more prominent trailing
component which is only $\sim 3\%$ linearly polarized. While the percentage of
circular polarization is $~{}18\%$ at 1.38 GHz (under the leading component),
it is almost negligible at 4.5 GHz. The PPA of the average profile follows an
S-shaped swing with a "notch" at the first peak (Figure 4) as also reported by
Kramer et al. (1999). This pulsar shows long–term profile variations (Kramer
et al., 1999), which have been attributed to improper polarization calibration
in the past (van Straten, 2013), even though more recent studies suggest that
intrinsic effects may likely be at play (Liu et al., 2015; Padmanabh et al.,
2021). By using the measured flux densities given in Table LABEL:tb:obs, we
measure a spectral index of $-1.43\pm 0.24$, consistent (within errors) with
the published value of $-1.7\pm 0.1$ (Kramer et al., 1998). Since the observed
flux density of this pulsar is known to vary significantly due to diffractive
scintillation (Liu et al., 2015; Shannon et al., 2014), flux density
measurements over multiple epochs are needed for accurate spectral index
measurements.

We do not see evidence for giant pulses. However, the data set consists of
very bright pulses of (S/N$\approx 30$), which are $\approx\times 5$ the mean
S/N. At 4.5 GHz, some bright single pulses appear at phase $\approx 0.53$
(phase $\approx 0.03$ away from the main component) as shown in Figure 1.

The brightest single pulses of PSR J1022+1001 at 1.38 GHz are highly linearly
polarized, as also noted in other studies (Liu et al., 2015). We find that the
single–pulse linear polarization distribution of PSR J1022+1001 at 1.38 GHz
shows two peaks (see Figure 4). Single pulses associated with the leading
component of the average profile tend to have a low linear polarization
fraction (average $L/I\approx 0.3$) compared to the single pulses associated
with the second peak (average $L/I\approx 0.7$). As evident from Table
LABEL:tb:pol and Figure 4, in general, the fraction of linear polarization is
reduced at higher frequencies. The brightest single pulses at 4.5 GHz (of $\rm
S/N>5$ which we consider in our analysis) show higher fractional $L/I$
compared to the average profile.

The rms timing residual due to pulse–to–pulse jitter at 1.38 GHz, scaled to
one hour, is $135\pm 29$ ns. This is consistent (within the 2-$\sigma$ error)
with the $67\pm 9$ reported by Feng et al. (2020) or $\approx 700$ ns for
1-min integration ($\approx 90$ ns for one hour) reported by Liu et al. (2015)
at 1.38 GHz, but lower than the Shannon et al. (2014) value of $290\pm 15$ ns
and $269\pm 4$ ns of Lam et al. (2019). The level of jitter at 4.5 GHz, scaled
to one hour, is 117$\pm$52 ns, which is consistent with that at 1.38 GHz
within errors.

Following Liu et al. (2015) and Feng et al. (2020), we measured the rms timing
residuals of single pulses from the leading and trailing component separately.
Considering all single pulses with $\rm S/N>5$ results in rms timing residuals
of $\approx 70\,\mu$s and $\approx 43\,\mu$s for leading and trailing
components respectively. For TOA calculation, the template profile was created
by averaging over all pulses of $\rm S/N>5$ within each component. Liu et al.
(2015) finds rms timing residuals of $4.3\,\mu$s for averages of 30 trailing
component sub-pulses of $\rm S/N>5$, which converts to
$4.3\times\sqrt{30}\approx 23\,\mu$s at the single–pulse level. The reason for
higher rms timing residuals in our dataset could be the higher levels of
jitter compared to Liu et al. (2015) possibly due to the pulsar showing
variable jitter noise or the difference in radiometer noise.

We also measured the rms timing residuals of single pulses within S/N ranges
of $5<\rm S/N<15$, $15<\rm S/N<25$ and S/N $>25$. As shown in Figure 7, the
rms timing residual using the brightest pulses ($\rm S/N>25$) is $\approx
22\,\mu$s. For TOA calculation, the template profile was created by averaging
over all pulses within each S/N range. Scaling the rms of $\approx 2.09\,\mu$s
from 1000-pulse averages (Table LABEL:tb:jitter) simply as
$\sigma_{obs}\propto 1/\sqrt{N_{p}}$ (e.g. Shannon & Cordes, 2012), where
$N_{p}$ is the number of pulses in the average profile, the expected
single–pulse RMS is $\approx 66\,\mu$s. So timing just the brightest pulses
gives a factor of $\approx 3$ improvement. The vast majority of pulses (99%)
within $\rm S/N>25$ range appear at the trailing component. Using single
pulses within the S/N ranges of $5<\rm S/N<15$ and $15<\rm S/N<25$ yields rms
timing residuals of 69 and 74 $\mu$s respectively. The higher rms residuals in
pulses within the $15<\rm S/N<25$ range could be due to a small percentage
(3.7%) of pulses occurring at the leading component while the rest occur at
the trailing component. When these trailing edge pulses are removed, the rms
timing residuals within $15<\rm S/N<25$ reduce to 28 $\mu$s.

###  4.2 PSR J1713+0747

PSR J1713+0747 is a binary pulsar which shows relativistic effects (Camilo et
al., 1994; Splaver et al., 2005). PSR J1713+0747 shows a multi–component
average profile with a significant linear polarization ($\approx 37\%$) at the
peak at 1.38 GHz, consistent with previous studies (e.g. Yan et al., 2011).
This is reduced but still significant ($\approx 30\%$) at 4.5 GHz. We measure
a spectral index of $-1.69\pm 0.38$ using the flux densities obtained at 1.38
GHz and 4.5 GHz (see Table LABEL:tb:obs), which is consistent (within errors)
with the published value of $-1.5\pm 0.1$ (Kramer et al., 1998). The more
prominent outer components at 4.5 GHz compared to at 1.38 GHz (at phases
$\approx$0.42 and 0.47), also reported in Kijak et al. (1997), are evident in
our data (see Figure 1). The faint emission component at phase $\approx$0.63
at 1.38 GHz is also more prominent at 4.5 GHz. Previous studies have reported
a jitter contribution of $\approx 35\pm 0.8$ ns to rms timing residuals
residuals in one hour of integration of this pulsar (Shannon et al., 2014; Liu
et al., 2016), consistent (within the $2\sigma$ error) with the $45\pm 7$ ns
that we find at 1.38 GHz. At 4.5 GHz, the lack of correlation (see the top two
panels of Figure 6) and the large error in the quadrature difference method
point to jitter being not constrained.

Although, the data set consists of very bright pulses of S/N$\approx 20$, we
do not see evidence for giant pulses. The fraction of single pulses detected
at 4.5 GHz is very small, which stops us from performing detailed analysis.

The PPA of the average profile shows 90∘ jumps (see Figure 5) that indicate
the presence of OPMs (Xilouris et al., 1998). Indeed OPMs have been detected
in this pulsar at 1.38 GHz in Large European Array for Pulsars (LEAP) data
(Liu et al., 2016). A faint offset in the PPA distribution of single pulses at
phase $\approx 0.5$ is noted in Figure 5, hinting at an OPM for this pulsar.

###  4.3 PSR B1855+09

PSR B1855+09 is a relativistic binary which has both main pulse (MP) and an
interpulse (IP), with both showing a double–peaked structure and weak
($\approx 7\%$ and $\approx 2\%$ in the MP and IP respectively) linear
polarization consistent with previous studies (Kramer et al., 1998). The PPA
of the average profile in the MP does not follow an S–shaped swing. The
discontinuities in the PPA curve of the IP is non–orthogonal and consistent
with (Xilouris et al., 1998) even though some studies have noted orthogonal
transitions (Yan et al., 2011). Based on the measured flux densities at 430
MHz and 1.38 GHz, we measure a spectral index of $-0.6\pm 0.1$, slightly
flatter than the published value of $-1.3\pm 0.2$ (Kramer et al., 1998). The
flux density measurement at 2 GHz is unusually low, likely due to the very low
S/N of the data set or scintillation, and therefore has not been included in
the spectral index measurement. We also note that our data are only
polarization calibrated, and therefore these flux density measurements may not
be precise. Therefore we again emphasise that multi–epoch flux density
measurements may give a better estimate of the the spectral index. Also note
that the MSP spectral index is not well constrained (Bates et al., 2013;
Kuniyoshi et al., 2015). We find that the rms timing residual due to jitter,
scaled to one hour, is $60\pm 25$ ns ($70\pm 11\,\mu$s at the single–pulse
level). This is roughly consistent within the $3\sigma$ confidence interval of
$90-115\,\mu$s published in Lam et al. (2019).

Single pulses of PSR B1855+09 at 1.38 GHz also show higher fractional $L/I$
compared to the average profile. We also note that even though the flux
density at 430 MHz is high, the expected S/N of a single pulse is low and
therefore single pulses are not detected at these low frequencies possibly due
to the high sky temperature. For Arecibo observations (parameters listed in
Table LABEL:tb:obs), the expected S/N, for B1855+09 at 430 MHz, calculated
from the radiometer equation,

$S/N=\frac{G\sqrt{N_{p}W\,\Delta\nu}\,S_{\rm peak}}{T_{\rm sys}}$ (6)

is $<5$. Here $S_{\rm peak}$ is the peak flux density (estimated using the
mean flux density and duty cycle) and we assume that the single-pulse width is
$W\approx 0.17$ ms, which is the mean single–pulse equivalent width of
B1855+09 at 1.38 GHz.

A few single pulses associated with the IP are detected with low S/N (example
shown in the bottom panel of Figure 1). However, due to their low S/N, only
single–pulses associated with the MP are used for detailed analysis.

##  5 Conclusion and discussion

We have studied the total intensity properties, single–pulse polarimetry and
pulse phase jitter of MSPs J1022+1001, J1713+0747, B1855+09 at multiple
frequencies. On the single–pulse front, we do not detect giant pulses, but we
still find bright single pulses at 1.38 GHz for all three MSPs. Overall,
single pulses with higher S/N seems to be narrower than the ones with smaller
S/N at 1.38 GHz, as shown in Figure 3.

As evident from Table LABEL:tb:pol, in general, the fraction of linear
polarization in the average profile is reduced at higher frequencies. This may
be consistent with similar findings of normal pulsars that show that average
profiles are depolarized at higher frequencies (e.g. Xilouris et al., 1996;
Johnston et al., 2006). This is because high frequency emission arises closer
to the pulsar surface where the magnetic field deviates from a dipole field
resulting in superposition of modes and hence significant depolarization (e.g.
McKinnon, 1997; Wang et al., 2015; van Straten & Tiburzi, 2017). Dai et al.
(2015) reported this trend for PSRs J1022+1001, J1713+0747 and B1855+09 using
1.38 GHz and 3 GHz observations. Our results confirm that this trend continues
at 4.5 GHz for PSRs J1022+1001 and J1713+0747. Single pulses from PSR
J1022+1001 at 4.5 GHz and PSR B1855+09 at 1.38 GHz show high fractional linear
polarization than the average profile. This could be because bright pulses
show more fractional polarization or due to depolarization from averaging.

In our data set the fraction of single pulses detected at 4.5 GHz is very
small, which prevents us from performing detailed analysis. Sensitive
observations are needed to better understand single–pulse phenomena at higher
frequencies. For example, the expected S/N at 4.5 GHz for a FAST–like
telescope (even though FAST cannot observe at frequencies > 2 GHz with current
system), with a gain of $\rm G=18\rm\,K/Jy$, $\rm T_{sys}=20$ K, $\rm BW=800$
MHz (Li & Pan, 2016), would be a factor of $\approx$ 6 better than our data.
Upcoming facilities with high sensitivity and frequency coverage should
facilitate multi–frequency single–pulse studies of pulsars.

The jitter levels ($\sigma_{j}$) we find are consistent for PSRs J1022+1002
(Liu et al., 2015; Feng et al., 2020; Parthasarathy et al., 2021), PSR
J1713+0747 (Shannon et al., 2014; Liu et al., 2016) and PSR B1855+09 (Lam et
al., 2019) at 1.38 GHz. But the fact that PSR J1022+1001 has lower values of
jitter recorded in other studies (Shannon et al., 2014; Lam et al., 2019)
could indicate variable jitter noise. The contribution to the pulsar noise
budget from jitter at 4.5 GHz and 1.38 GHz are consistent (within the
$2\sigma$ error) for PSR J1022+1001. The level of jitter from the cross-
correlation method ($\sigma_{\rm j,CCF}$) and the quadrature difference method
($\sigma_{j}$ in Table LABEL:tb:jitter) are comparable at 1.38 GHz. At 4.5 GHz
either the radiometer noise component is underestimated or jitter is
decorrelated across the band (e.g. Parthasarathy et al., 2021). Jitter noise
is expected to be lower at higher frequencies due to narrower profiles and
lower flux densities, though previous studies have shown that not all pulsars
follow this trend (Shannon et al., 2014; Lam et al., 2019). While some MSPs
used in PTAs show levels of jitter as low as tens of ns per hour, others
(which have large pulse widths and small periods) show levels of jitter
amounting to hundreds of ns per hour (Shannon et al., 2014; Parthasarathy et
al., 2021). Lam et al. (2019) has found that 43 pulsars out of a sample of 48
shows significant jitter and that 30 of those show significant frequency
dependence. The general notion is that jitter noise should be detectable in
all MSPs if observed with adequate sensitivity (Shannon et al., 2014), and
therefore accurately including jitter in noise modeling is important in high-
precision pulsar timing. It also represents a fundamental lower limit on the
timing precision achievable for individual MSPs.

We also find that timing only the brightest pulses of PSR J1022+1001 results
in an improvement in rms timing residuals by a factor of $\approx$3 than what
is expected by scaling timing properties derived from all single pulses.
However, given relative paucity of such bright pulses, the timing may not be
improved when compared to fully integrating all available data in a fixed
observation time. We note that previous studies have found that selective
timing does not yield a significant improvement in timing (Osłowski et al.,
2014; Feng et al., 2020). The first of these studies reported reduced
$\chi^{2}$ of the timing residuals when selectively integrating only weaker
pulses, as well as a minor improvement to timing when using brighter pulses.
The reason for our improvement in timing could be because we only include the
brightest pulses with a large trailing component and remove the pulses that
have a large leading component and thereby selecting only the brightest pulses
which appear within the narrow phase range of the trailing component. Our very
short data set prevents us from timing using integrated profiles of selected
bright pulses, but this is an important avenue for further exploration.
Therefore, further investigation is needed to see if selective timing yields
an improvement for integrated profiles formed with bright pulses within a
given component of the profile.

## Data availability

The data underlying this article will be shared on reasonable request to the
corresponding author.

## acknowledgement

The authors thank Hector Hernandez, Phil Perrillat and other Arecibo
Observatory staff for help with scheduling, observation support and data
quality checks. The Arecibo Observatory is a facility of the National Science
Foundation operated under cooperative agreement by the University of Central
Florida in alliance with Yang Enterprises, Inc. and Universidad Metropolitana.
The authors also thank the anonymous referee for the constructive remarks on
this manuscrpt. NP, BPP, and MAM are members of the NANOGrav Physics Frontiers
Center, supported by NSF award #1430284 and #2020265.

## References

  * Alam et al. (2021) Alam, MD F. et al. 2021, ApJS, 252, 4A
  * Arzoumanian et al. (2018) Arzoumanian, Z., Brazier, A., Burke-Spolaor, S., et al. 2018, ApJS, 235, 37. doi:10.3847/1538-4365/aab5b0
  * Backer (1970) Backer, D. C. 1970, Nature, 228, 1297
  * Backer (1973) Backer, D. C. 1973, ApJ, 182, 245
  * Backer et al. (1976) Backer, D. C., Rankin, J. M., & Campbell, D. B. 1976, Nat, 263, 202
  * Backer, & Rankin (1980) Backer, D. C., & Rankin, J. M. 1980, ApJS, 42, 143
  * Bates et al. (2013) Bates, S. D., Lorimer, D. R., & Verbiest, J. P. W. 2013, MNRAS, 431, 1352. doi:10.1093/mnras/stt257
  * Britton (2000) Britton, M. C. 2000, ApJ, 532, 1240
  * Brook et al. (2018) Brook, P. R., Karastergiou, A., McLaughlin, M. A., et al. 2018, ApJ, 868, 122. doi:10.3847/1538-4357/aae9e3
  * Cocke (1973) Cocke, W. J. 1973, ApJ, 184, 291. doi:10.1086/152326
  * Cordes & Shannon (2010) Cordes, J. M., & Shannon, R. M. 2010, arXiv:1010.3785
  * Cordes et al. (1978) Cordes, J. M., Rankin, J., & Backer, D. C. 1978, ApJ, 223, 961
  * Cordes et al. (1990) Cordes, J. M., Weisberg, J. M., & Hankins, T. H. 1990, AJ, 100, 1882. doi:10.1086/115644
  * Dai et al. (2015) Dai, S., Hobbs, G., Manchester, R. N., et al. 2015, MNRAS, 449, 3223. doi:10.1093/mnras/stv508
  * De et al. (2016) De, K., Gupta, Y., & Sharma, P. 2016, ApJ, 833, L10. doi:10.3847/2041-8213/833/1/L10
  * Drake & Craft (1968) Frake, F. D., & Craft, H. D. 1968, Nature, 220, 231
  * Desvignes et al. (2016) Desvignes, G. et al. 2016, MNRAS, 458, 3341
  * Downs & Reichley (1983) Downs, G. S. & Reichley, P. E. 1983, ApJS, 53, 169. doi:10.1086/190890
  * Edwards et al. (2006) Edwards, R. T., Hobbs, G. B., & Manchester, R. N. 2006, MNRAS, 372, 1549. doi:10.1111/j.1365-2966.2006.10870.x
  * Enoto et al. (2021) Enoto, T., Terasawa, T., Kisaka, S., et al. 2021, Science, 372, 187. doi:10.1126/science.abd4659
  * Everett & Weisberg (2001) Everett, J. E. & Weisberg, J. M. 2001, ApJ, 553, 341. doi:10.1086/320652
  * Feng et al. (2020) Feng, Y., Hobbs, G., Li, D., et al. 2020, arXiv:2012.06709
  * Gentile et al. (2018) Gentile, P. A., McLaughlin, M. A., Demorest, P. B., et al. 2018, ApJ, 862, 47. doi:10.3847/1538-4357/aac9c9
  * Gil & Snakowski (1990) Gil, J. A., & Snakowski, J. K. 1990, AAP, 234, 269
  * Gil et al. (2002) Gil, J. A., Melikidze, G. I., & Mitra, D. 2002, A&A, 388, 235. doi:10.1051/0004-6361:20020473
  * Goldreich & Julian (1969) Goldreich, P. & Julian, W. H. 1969, ApJ, 157, 869. doi:10.1086/150119
  * Camilo et al. (1994) Camilo, F., Foster, R. S., & Wolszczan, A. 1994, ApJ, 437, L39. doi:10.1086/187677
  * Hankins & Eilek (2007) Hankins, T. H. & Eilek, J. A. 2007, ApJ, 670, 693. doi:10.1086/522362
  * Haslam et al. (1982) Haslam, C. G. T., Salter, C. J., Stoffel, H., et al. 1982, A&AS, 47, 1
  * Heiles & Campbell (1970) Heiles, C. & Campbell, D. B. 1970, Nature, 226, 529. doi:10.1038/226529a0
  * Hobbs et al. (2006) Hobbs, G. B., Edwards, R. T., & Manchester, R. N. 2006, MNRAS, 369, 655. doi:10.1111/j.1365-2966.2006.10302.x
  * Hotan et al. (2004) Hotan, A. W., Bailes, M., & Ord, S. M. 2004, MNRAS, 355, 941. doi:10.1111/j.1365-2966.2004.08376.x
  * Hotan et al. (2006) Hotan, A. W., Bailes, M., & Ord, S. M. 2006, MNRAS, 369, 1502. doi:10.1111/j.1365-2966.2006.10394.x
  * Hotan et al. (2004) Hotan, A. W., van Straten, W., & Manchester, R. N. 2004, PASA, 21, 302
  * Jenet et al. (1998) Jenet, F. A., Anderson, S. B., Kaspi, V. M., Prince, T. A., & Unwin, S. C. 1998, ApJ, 498, 365
  * Jenet et al. (2001) Jenet, F. A., Anderson, S. B., & Prince, T. A. 2001, ApJ, 546, 394
  * Johnston et al. (2006) Johnston, S., Karastergiou, A., & Willett, K. 2006, MNRAS, 369, 1916. doi:10.1111/j.1365-2966.2006.10440.x
  * Kalapotharakos et al. (2021) Kalapotharakos, C., Wadiasingh, Z., Harding, A. K., et al. 2021, ApJ, 907, 63. doi:10.3847/1538-4357/abcec0
  * Karuppusamy et al. (2012) Karuppusamy, R., Stappers, B. W., & Lee, K. J. 2012, A&A, 538, A7. doi:10.1051/0004-6361/201117667
  * Kerr et al. (2020) Kerr, M. et al. 2020, PASA, 37, 20K
  * Kijak et al. (1997) Kijak, J., Kramer, M., Wielebinski, R., et al. 1997, A&A, 318, L63
  * Knight (2007) Knight, H. S. 2007, MNRAS, 378, 723. doi:10.1111/j.1365-2966.2007.11810.x
  * Kramer et al. (1998) Kramer, M., Xilouris, K. M., Lorimer, D. R., et al. 1998, ApJ, 501, 270
  * Kramer et al. (1999) Kramer, M., Xilouris, K. M., Camilo, F., et al. 1999, ApJ, 520, 324. doi:10.1086/307449
  * Kuniyoshi et al. (2015) Kuniyoshi, M., Verbiest, J. P. W., Lee, K. J., et al. 2015, MNRAS, 453, 828. doi:10.1093/mnras/stv1604
  * Lam et al. (2019) Lam, M. T., McLaughlin, M. A., Arzoumanian, Z., et al. 2019, ApJ, 872, 193. doi:10.3847/1538-4357/ab01cd
  * Liu et al. (2011) Liu, K., Verbiest, J. P. W., Kramer, M., et al. 2011, MNRAS, 417, 2916. doi:10.1111/j.1365-2966.2011.19452.x
  * Li & Pan (2016) Li, D. & Pan, Z. 2016, Radio Science, 51, 1060
  * Liu et al. (2015) Liu, K., Karuppusamy, R., Lee, K. J., et al. 2015, MNRAS, 449, 1158. doi:10.1093/mnras/stv397
  * Liu et al. (2016) Liu, K., Bassa, C. G., Janssen, G. H., et al. 2016, MNRAS, 463, 3239. doi:10.1093/mnras/stw2223
  * Lorimer & Kramer (2012) Lorimer, D. R. & Kramer, M. 2012, Handbook of Pulsar Astronomy, by D. R. Lorimer , M. Kramer, Cambridge, UK: Cambridge University Press, 2012
  * Lyne & Smith (1989) Lyne, A. G. & Smith, F. G. 1989, MNRAS, 237, 533. doi:10.1093/mnras/237.3.533
  * Lyne et al. (2010) Lyne, A., Hobbs, G., Kramer, M., et al. 2010, Science, 329, 408. doi:10.1126/science.1186683
  * Machabeli & Usov (1979) Machabeli, G. Z. & Usov, V. V. 1979, Soviet Astronomy Letters, 5, 238
  * Maron et al. (2004) Maron, O., Kijak, J., & Wielebinski, R. 2004, A&A, 413, L19. doi:10.1051/0004-6361:20031707
  * McKee et al. (2019) McKee, J. W., Stappers, B. W., Bassa, C. G., et al. 2019, MNRAS, 483, 4784. doi:10.1093/mnras/sty3058
  * McKinnon (1997) McKinnon, M. M. 1997, ApJ, 475, 763. doi:10.1086/303542
  * McKinnon & Stinebring (2000) McKinnon, M. M. & Stinebring, D. R. 2000, ApJ, 529, 435. doi:10.1086/308264
  * Melrose (1992) Melrose, D. B. 1992, The Magnetospheric Structure and Emission Mechanisms ofRadio Pulsars, IAU Colloquium 128, 306
  * Melrose (1995) Melrose, D. B. 1995, Journal of Astrophysics and Astronomy, 16, 137
  * Melrose & Gedalin (1999) Melrose, D. B. & Gedalin, M. E. 1999, ApJ, 521, 351. doi:10.1086/307539
  * Melrose et al. (2021) Melrose, D. B., Rafat, M. Z., Mastrano, A. 2021, MNRAS, 500, 4530
  * Mitra et al. (2009) Mitra, D., Gil, J., & Melikidze, G. I. 2009, ApJ, 696, L141. doi:10.1088/0004-637X/696/2/L141
  * Mitra & Rankin (2011) Mitra, D. & Rankin, J. M. 2011, ApJ, 727, 92
  * Noutsos et al. (2015) Noutsos, A., Sobey, C., Kondratiev, V. I., et al. 2015, A&A, 576, A62. doi:10.1051/0004-6361/201425186
  * Osłowski et al. (2013) Osłowski, S., van Straten, W., Demorest, P., & Bailes, M. 2013, MNRAS, 430, 416
  * Osłowski et al. (2014) Osłowski, S., van Straten, W., Bailes, M., et al. 2014, MNRAS, 441, 3148. doi:10.1093/mnras/stu804
  * Padmanabh et al. (2021) Padmanabh, P. V., Barr, E. D., Champion, D. J., et al. 2021, MNRAS, 500, 1178. doi:10.1093/mnras/staa3174
  * Parthasarathy et al. (2021) Parthasarathy, A., Bailes, M., Shannon, R. M., et al. 2021, MNRAS, 502, 407. doi:10.1093/mnras/stab037
  * Perera et al. (2018) Perera, B. B. P. et al. 2018, MNRAS, 478, 218
  * Perera et al. (2019) Perera, B. B. P. et al. 2019, MNRAS, 490, 4666
  * Radhakrishnan & Cooke (1969) Radhakrishnan, V. & Cooke, D. J. 1969, Astrophys. Lett., 3, 225
  * Rankin (1983a) Rankin, J. M. 1983, ApJ, 274, 333
  * Rankin (1983b) Rankin, J. M. 1983, ApJ, 274, 359
  * Rankin (1986) Rankin, J. M. 1986, ApJ, 301, 901. doi:10.1086/163955
  * Romani & Johnston (2001) Romani, R. W. & Johnston, S. 2001, ApJ, 557, L93. doi:10.1086/323415
  * Shannon et al. (2014) Shannon, R. M., Osłowski, S., Dai, S., et al. 2014, MNRAS, 443, 1463. doi:10.1093/mnras/stu1213
  * Splaver et al. (2005) Splaver, E. M., Nice, D. J., Stairs, I. H., et al. 2005, ApJ, 620, 405. doi:10.1086/426804
  * Reardon et al. (2016) Reardon, D. J., Hobbs, G., Coles, W., et al. 2016, MNRAS, 455, 1751. doi:10.1093/mnras/stv2395
  * Ruderman & Sutherland (1975) Ruderman, M. A. & Sutherland, P. G. 1975, ApJ, 196, 51. doi:10.1086/153393
  * Shannon & Cordes (2012) Shannon, R. M., & Cordes, J. M. 2012, arXiv:1210.7021
  * Simard-Normandin & Kronberg (1980) Simard-Normandin, M. & Kronberg, P. P. 1980, ApJ, 242, 74. doi:10.1086/158445
  * van Straten (2004) van Straten, W. 2004, ApJS, 152, 129
  * van Straten (2006) van Straten, W. 2006, ApJ, 642, 1004
  * van Straten & Bailes (2011) van Straten, W., & Bailes, M. 2011, PASA, 28, 1
  * van Straten (2013) van Straten, W. 2013, ApJS, 204, 13. doi:10.1088/0067-0049/204/1/13
  * van Straten & Tiburzi (2017) van Straten, W., & Tiburzi, C. 2017, ApJ, 835, 293. doi:10.3847/1538-4357/835/2/293
  * Wahl et al. (2021) Wahl, H. M., McLaughlin, M., Gentile, P. A., et al. 2021, arXiv:2104.05723
  * Wang et al. (2015) Wang, P. F., Wang, C., & Han, J. L. 2015, MNRAS, 448, 771. doi:10.1093/mnras/stu2765
  * Yan et al. (2011) Yan, W. M., Manchester, R. N., van Straten, W., et al. 2011, MNRAS, 414, 2087. doi:10.1111/j.1365-2966.2011.18522.x
  * You et al. (2012) You, X. P., Coles, W. A., Hobbs, G. B., et al. 2012, MNRAS, 422, 1160. doi:10.1111/j.1365-2966.2012.20688.x
  * Xilouris et al. (1996) Xilouris, K. M., Kramer, M., Jessner, A., et al. 1996, A&A, 309, 481
  * Xilouris et al. (1998) Xilouris, K. M., Kramer, M., Jessner, A., et al. 1998, ApJ, 501, 286. doi:10.1086/305791
  * Xu et al. (2021) Xu, H., Huang, Y. X., Burgay, M., et al. 2021, The Astronomer’s Telegram, 14642
  * Zhelezniakov & Shaposhnikov (1979) Zhelezniakov, V. V. & Shaposhnikov, V. E. 1979, Australian Journal of Physics, 32, 49. doi:10.1071/PH790049
  * Zhu et al. (2015) Zhu, W. W., Stairs, I. H., Demorest, P. B., et al. 2015, ApJ, 809, 41. doi:10.1088/0004-637X/809/1/41


# Explainable deepfake and spoofing detection:
an attack analysis using SHapley Additive exPlanations

###### Abstract

Despite several years of research in deepfake and spoofing detection for
automatic speaker verification, little is known about the artefacts that
classifiers use to distinguish between bona fide and spoofed utterances. An
understanding of these is crucial to the design of trustworthy, explainable
solutions. In this paper we report an extension of our previous work to better
understand classifier behaviour to the use of SHapley Additive exPlanations
(SHAP) to attack analysis. Our goal is to identify the artefacts that
characterise utterances generated by different attacks algorithms. Using a
pair of classifiers which operate either upon raw waveforms or magnitude
spectrograms, we show that visualisations of SHAP results can be used to
identify attack-specific artefacts and the differences and consistencies
between synthetic speech and converted voice spoofing attacks.

##  1 Introduction

Deep neural networks (DNNs) are now the de facto approach to spoofing
detection for automatic speaker verification [1, 2, 3, 4]. Spoofing
countermeasures usually operate upon hand-crafted spectro-temporal features
extracted from the input utterance [5, 6, 7, 8, 9, 10]. Lately, so-called raw
end-to-end (E2E) solutions [11, 12, 13, 14, 15] are beginning to outperform
the more traditional approaches. By operating directly upon raw audio
waveforms, such systems have greater potential to capture the tell-tale signs
of spoofing attacks, e.g. speech synthesis, converted voice and replay, the
artefacts of which might not be captured in handcrafted features that are
often designed for other tasks.

Whatever the approach, there is usually an interest to understand _how_ or
_why_ a given classifier arrives at the scores or the decisions it produces.
Such an understanding is crucial to the design of trustworthy systems and key
to the trend toward explainable artificial intelligence (xAI) [16, 17]. So
far, while we know that different spoofing detection solutions can rely upon
different temporal or spectral intervals [18, 11, 19], we know surprisingly
little else, despite several years of research.

An explanation of the cues used by a classifier in distinguishing between bona
fide and spoofed speech is of interest not just in terms of curiosity.
Explanations can be crucial to some particular use case scenarios, e.g. those
involving forensics, but they might also reveal opportunities to improve
performance through feature engineering and/or classifier design. Explanations
might also be particularly relevant in the case of now-dominant deep neural
networks, increasingly opaque, black-box spoofing detection solutions with no
obvious means to help explain classifier behaviour, scores or decisions.

Motivated by related explainability studies in other speech tasks [20, 21,
22], we have investigated the use of SHapley Additive exPlanations (SHAP) [23]
to explore and compare the behaviour of DNN-based solutions to spoofing
detection [24]. SHAP is a feature attribution approach to explainability. A
so-called SHAP value is calculated and assigned individually to each and every
input feature. For image-related tasks, the feature can be an image pixel. For
speech-related tasks, the feature can be a spectro-temporal magnitude estimate
or, as we show later in this paper, a sample of the input time-domain
waveform. SHAP values reflect the difference in classifier output when derived
either with or without the use of each feature. The SHAP value hence indicates
the influence or importance of each input feature upon the classifier output.

SHAP values can be readily visualised with heat maps of the same dimension as
the input image or spectro-temporal decomposition. We followed this approach
in our previous work to explore differences in _classifier behaviour_ and
confirmed that different solutions exploit different temporal intervals and
spectral sub-bands to distinguish between bona fide and spoofed utterances,
behaviour that cannot be derived from classifier scores or decisions alone.
The current article reports our work to extend the study to _attack analysis_.

Here we are interested to learn more about the attacks themselves, namely the
specific artefacts that characterise each attack or class of attacks, e.g.,
algorithms such as text-to-speech (TTS) or voice conversion (VC). We focus
upon the training data in the ASVspoof 2019 Logical Access (LA) database and
use a pair of classifiers trained using spoofed data generated with specific
attack algorithms to identify the characteristic artefacts. We found that this
approach is critical to the identification of consistent artefacts in training
data for which ground truth attack algorithms are available. This same
approach, though, is obviously impracticable in the case of evaluation data,
or spoofing detection in the wild, where the ground truth label of the
spoofing attack algorithm is, by default, unknown. In this case,
explainability proves to be considerably more challenging. Nonetheless, there
are similarities between the characteristics of different known attacks,
similarities which provide a starting point for explainability in the wild.

The new contributions in this work are: (i) the first application of SHAP
analysis to an end-to-end spoofing detection solution which operates upon raw
waveform inputs; (ii) an attack analysis in terms of explainability using two
different classifiers, one that operates upon spectro-temporal magnitude
estimates, the other upon raw waveform inputs; (iii) explainability
visualisations for both classifier inputs.

##  2 SHapley Additive exPlanations (SHAP)

We provide a brief description of SHAP analysis. In extending our previous
work [24], we also show how it can be applied to raw, temporal inputs and
provide example visualisations applied to classifiers that operate upon either
spectro-temporal decompositions or raw waveforms.

###  2.1 Definition

SHAP [23] is used to explain the contribution of each individual input feature
to the output of a given model $f(x)$, here a spoofing detection classifier.
The contribution of each feature is reflected via SHAP values $\phi$, a
reflection of the difference in the model output when it is learned with or
without the inclusion of this particular feature. When $f(\mathbf{x})$ is a
complex model which takes considerable time to retrain, e.g. deep neural
network solutions, some simplifications are necessary to mitigate the need for
repetitive retraining. First, the input $\mathbf{x}$ is approximated by a
simplified feature
$\mathbf{x}^{\prime}=\\{x^{\prime}_{1},...,x^{\prime}_{D}\\}$, where
$x^{\prime}_{n}\in\\{0,1\\}$ implies either the absence ($0$) or presence
($1$) of the corresponding feature in $\mathbf{x}$, and $D$ is the feature
dimension. The original model $f(\mathbf{x})$ is then approximated with an
explanation model $g(\mathbf{x}^{\prime})$:

$f(\mathbf{x})\approx
g(\mathbf{x}^{\prime})=\phi_{0}+\sum_{n=1}^{D}\phi_{n}x^{\prime}_{n}$ (1)

where $\phi_{n}$ is the SHAP value for feature $x_{n}$,
$\phi_{0}=f(h_{x}(\mathbf{0}))$ and $h_{x}$ is a mapping function which
converts simplified features $\mathbf{x}^{\prime}$ to original features
$\mathbf{x}$, i.e., $\mathbf{x}=h(\mathbf{x}^{\prime})$. The explanation model
$g(\mathbf{x}^{\prime})$ is trained to approximate the original model output
$f(\mathbf{x})$ using the sum of SHAP values corresponding to the features for
which $x^{\prime}_{n}=1$ and coefficients $\phi_{n}$ in place of true SHAP
values. SHAP values $\phi_{n}$ are derived for each class and can be positive
and negative valued, indicating support (or not) for the corresponding class
hypothesis. Further information and details are available in [23, 25].

(a) waveform

(b) SHAP values for bona fide (green) and spoofed (red)

Figure 1: Time waveform and SHAP values for the `LA_T_2909480’ A03 text-to-
speech spoofed utterance `_Well, Scotland had better grow up, fast_ ' from the
ASVspoof 2019 LA database and a 1D-Res-TSSDNet spoofing detection model.

###  2.2 Example visualisations

Our first use of SHAP analysis applied to explainable spoofing detection [24]
was performed using spectro-temporal decompositions in the form of magnitude
spectrograms. This choice was motivated by their use to study easily the
influence of different time intervals and frequency sub-bands upon the outputs
of different classifiers, with the emphasis being upon characterising
_classifier behaviour_.

The goal of the current work relates more to _attack analysis_. Attack
artefacts, though, are undeniably the result of classifier behaviour; and they
are inseparable within the current framework. Accordingly, in a modest attempt
to decouple or marginalise classifier influences, the work reported in this
paper was also performed with two different classifiers.

Nowadays, a growing number of competitive systems operate directly upon raw
waveforms [26, 27, 28]. Since these approaches have potential to discover
artefacts that might be lost through hand-crafted pre-processing or feature
extraction, there is hence an interest to apply SHAP analysis to classifiers
operating upon raw waveforms, in addition to spectro-temporal decomposition.
Our expectation is that the use of two such different classifiers may help to
provide a more generic account of spoofing artefacts that is less specific to
the behaviour or limitations of a single feature representation and classifier
alone.

####  2.2.1 Raw waveforms

The waveform of a spoofed utterance $x[n]$ selected from the A03 spoofed
utterances in the ASVspoof 2019 LA training partition is depicted in Fig.
1(a). Using an arbitrary classifier which operates on the raw waveform, SHAP
values $\phi_{n}$ are obtained using Eq. 1 for each sample $x[n]$ and each
class: bona fide and spoofed.

We have found that SHAP values for each class are approximately negatively
symmetric, i.e. positive SHAP values for the bona fide class are correlated
with negative values for the spoof class, and vice versa. To avoid redundancy,
everywhere in this paper, we consider only positive SHAP values $\phi_{n}$
that favour either class. They are plotted in Fig. 1(b) for the same
utterance. Those for the bona fide class are illustrated in green whereas
those for the spoof class are inverted and plotted below the abscissa in red.

The plot shows that the classifier uses information distributed across the
full utterance, though SHAP values are mostly larger for speech segments and
lower for non-speech segments. The largest SHAP values are observed for speech
segments around 0.9 seconds and 1.5 seconds which are most supportive of
spoofed class hypotheses. Hence, while there is information everywhere in the
utterance, the classifier output is dominated by cues or artefacts located in
a specific temporal intervals.

####  2.2.2 Magnitude spectrograms

The spectrogram of the same utterance as in Fig. 1(a) is illustrated in Fig.
2(a). In similar fashion to [22, 24, 20], SHAP values $\phi_{(m,n)}$ for an
arbitrary classifier which operates upon such spectro-temporal magnitude
estimates $\mathbf{X}(m,n)$, where $m$ is the spectral bin and $n$ is the
frame index, are illustrated in Fig. 2(b). Once again, we illustrate positive
SHAP values only, which are now encoded in the colour intensity, still green
for the bona fide class and red for the spoof class.

As for the temporal plots in Fig. 1, we find that information is scattered
across the full spectrogram, though SHAP values are again greater for speech
intervals than for non-speech intervals. Since the SHAP values are now encoded
in the colour intensity, it is more difficult to identify from the
visualisation the specific temporal or spectral intervals that have the
greatest influence upon the classifier output, even if the most influential
temporal segments appear to be different to those highlighted in Fig. 1.

(a) magnitude spectrogram

(b) SHAP values

Figure 2: As for Fig. 1 except for a magnitude spectrogram input
representation and a 2D-Res-TSSDNet spoofing detection model. SHAP values in
(b) are dilated to improve their visualisation.

##  3 Experimental setup

In this section we describe the classifiers used for the remainder of this
work presented in this paper, specific implementation details and the pruning
of SHAP values to better highlight the temporal and/or spectral intervals that
bear the greatest influence upon the classifier output.

###  3.1 Spoofing detection classifiers

We used a pair of different classifiers, namely the 1D- and 2D-Res-TSSDNet
systems proposed in [13]. The motivation for using these two particular
systems are: they are freely available as open
source111https://github.com/ghuawhu/end-to-end-synthetic-speech-detection and
can hence be used to reproduce our results; they share a similar network
architecture consisting of several convolutional blocks with residual
connections, followed by a global max pooling layer and three fully-connected
layers, and both produce scores for bona fide and spoofed classes; they are
efficient and can be trained quickly; their relatively uncomplicated
structures should offer better potential for explainability.

The 1D-Res-TSSDNet has only 0.35M trainable parameters whereas the 2D-Res-
TSSDNet still has only 0.97M; our own solutions have a considerably greater
number [12, 11]. The core difference between the two classifiers is that the
1D-Res-TSSDNet is based upon 1D convolutional operations applied initially to
raw audio waveforms, whereas the 2D-Res-TSSDNet is based upon 2D convolutional
operations applied to magnitude spectrograms extracted using 20 ms Hamming
windows with a 10 ms shift and a 320-point FFT.

###  3.2 Implementation details

All experiments reported in this paper were conducted using the ASVspoof 2019
Logical Access (LA) database [29]. The training partition is used to update
network parameters whereas the evaluation partition is used as validation to
select the best model. System performance is estimated using the equal error
rate (EER). Since the goal of this work is to study spoofing artefacts and not
to learn the best performing classifier, we make no use of the development
partition for any experiments and use EER estimates only for the purposes of
model selection. Both systems are optimised by minimising the weighted cross-
entropy loss between spoofed and bona fide classes. Models are trained using
the Adam optimiser [30] with a learning rate of 0.001 and an exponential
learning rate decay of 0.95.

Instead of using fixed-length inputs resulting from concatenation or
truncation in order to create fixed-size mini-batches, both models operate
upon variable, whole-length utterances. This choice was found to be beneficial
in terms of explainability since it allows us to study more easily and
consistently the influence upon classifier ouput of information available in
specific temporal or spectral intervals. We found that concatenation in
particular can result in a classifier making different use of the same
information in repeated short intervals or sub-bands [24], a behaviour which
can complicate studies in explainability.

Variable-length input data is treated as described in [31] whereby training
utterances of similar length are zero-padded to give uniformly-sized mini-
batches. For experiments involving SHAP analysis, as opposed to model
updating, utterances are nonetheless fed to the classifiers without any such
padding. Though the use of variable-length inputs does degrade performance, we
stress again that this work targets explainability, not performance. For
reference, the EERs of the 1D-Res-TSSDNet and 2D-Res-TSSDNet are 6.90% and
4.28% respectively.

###  3.3 SHAP analysis and pruning

SHAP values222https://captum.ai/api/gradient_shap.html for a given input
feature are computed from background features (also called baselines) of
identically-sized zero-valued vectors. The number of randomly generated
examples per sample is set to 20. All experiments reported in the paper are
performed with the same random seed on a single NVIDIA GeForce RTX 3090 GPU.
Results are reproducible with the same random seed and GPU environment using
the same SHAP implementation and our scripts available online.333Our scripts
will be made available upon publication.

While SHAP values are derived for every time sample $n$ or every specto-
temporal sample $[m,n]$, we observed particularly high SHAP values for
typically small numbers of specific samples or intervals. Histograms of SHAP
values for both 1D and 2D classifiers depicted in Fig. 3 show approximate
negative exponential distributions and that the highest SHAP values are
accounted for by approximately 0.2% of samples. By selectively pruning the
input, we confirmed through experiments not reported here that samples with
these highest SHAP values have a dominant influence upon the classifier output
whereas samples with near-to-zero SHAP values have comparatively less
influence. Accordingly, all SHAP visualisations reported in the remainder of
this paper correspond to the 0.2% of samples with the highest, dominant SHAP
values. In similar fashion to [20], SHAP values in support of the spoofed
class (red only) are encoded in colour intensity superimposed upon waveform or
magnitude spectrogram plots.

(a) waveform feature

(b) spectrogram feature

Figure 3: Histograms of SHAP values for both classifiers and the utterance
shown in Fig. 1.

##  4 Attack analysis

We present an attack analysis which focuses upon spoofing attacks contained in
the ASVspoof 2019 LA training partition. It includes spoofed utterances
generated with 6 different algorithms including 4 TTS attacks (A01-A04) and 2
VC attacks (A05-A06).

Table 1:  Artefact description of attacks in ASVspoof 2019 LA train partition. | Found artefacts
---|---
Attack |  Algorithm |  Waveform |  Spectrogram
A01 |  TTS |  Vowels |  Lower frequency bands, leading 0.5s
A02 |  TTS |  Single dominant vowel |  Lower & higher frequency bands, unvoiced \s\
A03 |  TTS |  Less densely distributed in vowels |  Lower frequency bands
A04 |  TTS |  Non-speech, low energy speech segments (voice onsets and offsets) |  Full spectrum, unvoiced speech, clicks
A05 |  VC |  Voice onset, vowels |  Full spectrum, higher energy formant frequencies
A06 |  VC |  Speech distortion |  Lower frequency bands

(a) A01 attack utterance LA_T_3566209

(b) A02 attack utterance LA_T_1590397

(c) A03 attack utterance LA_T_2909480

(d) A04 attack utterance LA_T_5116902

(e) A05 attack utterance LA_T_3134909

(f) A06 attack utterance LA_T_5300749

Figure 4: SHAP values for the A01-A04 utterance `_Well, Scotland had better
grow up, fast_ ' and A05 and A06 utterance `_It raises a serious question
mark_ '. SHAP values superimposed upon either waveforms or spectrograms and
dilated in case of the latter, in similar fashion to Fig. 2.

For all reported results, the classifiers are trained using bona fide data and
_matched_ spoofed data only – SHAP analysis for an A01 attack utterance is
performed using classifiers trained only with A01 attack utterances, not the
full set of A01-A06 training attack utterances. We found that this approach
helps in the identification of consistent artefacts for each attack. Our aim
is to isolate or marginalise the influence of classifier behaviour to the
extent that this is possible. We accept that, within our current framework, it
is impossible to _fully_ isolate classifier behaviour from attack artefacts,
since the latter are always learned using a specific classifier. By learning
the classifier with matched attack utterances, we hope at least to observe
more clearly and consistently the artefacts that are specific to the given
attack, and not those learned from other attacks. This approach is consistent
with our goal of _attack analysis_.

The following analyses were preformed upon a random selection of 100
utterances for each training attack. Example results are illustrated in Fig. 4
which shows both waveforms and spectrograms with the corresponding
superimposed 0.2% highest intensity-encoded SHAP values. The illustrated
examples are specifically chosen to use the same utterance for TTS attacks and
the same utterance for VC attacks444Text inputs used in generating TTS attack
utterances are not necessarily consistent with the text content of VC attacks
[32].. A discussion of results for each attack is presented in the following
whereas a summary of the principal, consistent artefacts observed in each case
is presented in Table 1. These were identified from SHAP analysis and casual
listening test.

A01 is a neural network (NN) based TTS attack with a WaveNet [33] vocoder. A
waveform and spectrogram with superimposed SHAP values are illustrated in Fig.
4(a). We observed differences for 1D and 2D classifiers. For the 1D classifier
we found most artefacts in vowel segments though we could not identify a
particular vowel for which SHAP values are consistently the highest. Most
artefacts are found at low frequency bands. For a substantial number of
utterances, the 2D classifier identifies artefacts in the leading 0.5 seconds
of speech. This might be the result of A01 attacks having a consistently
shorter leading silence interval compared to bona fide utterances [34].

A02 is also an NN-based TTS attack, but with a WORLD [35] vocoder. Results are
illustrated in Fig. 4(b). Like the A01 attack, the 1D classifier finds
artefacts in vowel segments, such as the \o\ vowel in the given example. The
2D classifier finds artefacts mostly at lower frequencies and also at higher
frequency bands above 6 kHz. For A02 attacks and the 2D classifier, consistent
artefacts are identified for the unvoiced sound \s\segments.

A03 is a different NN-based TTS attack also with a WORLD vocder. Results are
shown in Fig. 4(c). For the 1D classifier, artefacts are found mostly in vowel
segments, but are less densely distributed compared to A01 and A02 (relatively
fewer dark-red points). The reason might be that artefacts in A03 attacks are
located in particular samples which are different with the neighborhood ones.
While we found that artefacts lie mostly at lower frequencies in the case of
the 2D classifier, we did not succeed in identifying artefacts within
consistent speech segments.

A04 is a waveform concatenation TTS attack. Results are shown in Fig. 4(d).
For some A04 attacks we found artefacts to lie within leading non-speech and
low energy speech segments (onsets and offsets of speech), a characteristic
which differentiates A04 from the other TTS attacks. These observations may
correspond to the use of waveform concatenation and may explain why such
attacks generated with the MaryTTS platform [36] were initially challenging to
detect [37]. The 2D classifier uses cues throughout the full spectrogram.
Other consistent artefacts were found in unvoiced segments and click sounds
(around 1.2 and 3.5 seconds for the given example).

A05 is an NN-based VC attack. Results are shown in Fig. 4(e). In addition to
more dominant artefacts in vowels segments (around 2.0 seconds in the given
example), the 1D classifier also finds consistent artefacts in lower energy
voice onset segments (around 1.6 seconds). Similar to the TTS A04 attack, the
2D classifier finds artefacts across the full spectrum, rather than specific
sub-bands. Nonetheless, higher SHAP values correspond generally to lower
frequencies and higher energy formant frequencies (around 2 and 4 kHz).

A06 is transfer-function-based VC attack [38]. Results are shown in Fig. 4(f).
For the 1D classifier, we observed temporal intervals with high SHAP values to
correspond to noticeably distorted, speech sounds. These seem to correspond to
variations in unnaturally high velocity. Using the 2D classifier, artefacts
are found mostly at lower frequencies below 3 kHz. There were not found to
correspond consistently to any particular speech sounds other than high-energy
segments and neither do they correspond consistently to the distortions
identified using the 1D classifier.

##  5 Discussion and conclusion

Reported in the paper is our use of SHAP analysis to characterise the six
synthetic speech and converted voice spoofing attacks that make up the
training partition of the ASVspoof 2019 LA database. The analysis, performed
with a pair of classifiers operating upon either time-domain waveforms or
spectro-temporal magnitude estimates show the importance of analysing
artefacts using different representations and classifiers; they identify
different artefacts. While the artefacts are also attack-specific, we
nonetheless observed consistencies in attack classes and algorithms.

Three of the four synthetic speech attacks exhibit temporally-distinct
artefacts located within vowel segments. These three attacks are all neural
network based text-to-speech algorithms. The fourth synthetic speech attack,
by contrast a waveform concatenation based approach, exhibits different
artefacts located within lower energy segments. Differences in artefacts are
possibly indicative of the spoofing attack algorithm. Differences were also
observed for the two voice conversion algorithms, with one showing artefacts
during voice onset segments and distributed throughout the spectrum, and the
other corresponding to notable distortions and lower frequencies. There are,
however, some consistencies between the artefacts produced by both converted
voice and synthetic speech attacks, particularly artefacts in low frequency
bands. While differences in artefacts may paint an ominous picture for the
designers of spoofing detection systems, consistencies are encouraging since
they offer some insights as to why, as shown in broader literature,
classifiers trained using spoofed utterances generated using only a small set
of spoofing algorithms are nonetheless effective in the face of utterances
generated with different spoofing algorithms.

Limitations of the current work relate to the use of classifiers trained using
matching spoofed utterances. This luxury is unavailable for practical use
cases for which the nature of spoofing attacks is unknown and unpredictable.
Explainability for _ideal_ scenarios involving _known_ attacks is already
extremely challenging. Even in these conditions the identification of
consistent trends is difficult. We are now working to apply the same technique
to different, _unknown_ spoofing attacks in a manner that nonetheless
identifies consistent artefacts. This work, using utterances contained in the
ASVspoof 2019 LA evaluation partition generated by 13 different attacks, is
proving even more challenging.

One current direction for our future work involves the identification of
attack-specific artefacts using matched training data and then their re-
identification using the same attack-specific classifiers used in this work.
Explainable spoofing detection solutions could then make use of SHAP
visualisations as presented in this paper, but also comparisons to artefacts
previously identified for known attack algorithms. Another direction for
future work involves the linking of artefacts to specific attack algorithm
components, e.g. a particular representations, neural architectures or
vocoders.

##  6 Acknowledgements

This work is supported by the TReSPAsS-ETN project funded from the European
Union’s Horizon 2020 research and innovation programme under the Marie
Skłodowska-Curie grant agreement No.860813. It is also supported by the
ExTENSoR project funded by the French Agence Nationale de la Recherche (ANR).

## References

  * [1] G. Lavrentyeva, S. Novoselov, A. Tseren, M. Volkova, A. Gorlanov and A. Kozlov,  ``STC antispoofing systems for the ASVspoof2019 challenge,''  in Proc. Interspeech, 2019, pp. 1033–1037.
  * [2] V. Peddinti, D. Povey and S. Khudanpur,  ``A time delay neural network architecture for efficient modeling of long temporal contexts,''  in Proc. Interspeech, 2015, pp. 3214–3218.
  * [3] C.-I. Lai, N. Chen, J. Villalba and N. Dehak,  ``ASSERT: Anti-spoofing with squeeze-excitation and residual networks,''  in Proc. Interspeech, 2019, pp. 1013–1017.
  * [4] D. S. Park, W. Chan, Y. Zhang, C.-C. Chiu, B. Zoph, E. D. Cubuk et al.,  ``SpecAugment: A simple data augmentation method for automatic speech recognition,''  in Proc. Interspeech, 2019, pp. 2613–2617.
  * [5] A. Nautsch, X. Wang, N. Evans, T. Kinnunen, V. Vestman, M. Todisco et al.,  ``ASVspoof 2019: spoofing countermeasures for the detection of synthesized, converted and replayed speech,''  IEEE Transactions on Biometrics, Behavior, and Identity Science, vol. 3, no. 2, pp. 252–265, 2021.
  * [6] X. Li, N. Li, C. Weng, X. Liu, D. Su, D. Yu et al.,  ``Replay and synthetic speech detection with Res2Net architecture,''  in Proc. ICASSP, 2021, pp. 6354–6358.
  * [7] Y. Zhang, F. Jiang and Z. Duan,  ``One-class learning towards synthetic voice spoofing detection,''  IEEE Signal Processing Letters, vol. 28, pp. 937–941, 2021.
  * [8] T. Chen, A. Kumar, P. Nagarsheth, G. Sivaraman and E. Khoury,  ``Generalization of audio deepfake detection,''  in Proc. Odyssey, 2020, pp. 132–137.
  * [9] A. Tomilov, A. Svishchev, M. Volkova, A. Chirkovskiy, A. Kondratev and G. Lavrentyeva,  ``STC Antispoofing Systems for the ASVspoof2021 Challenge,''  in Proc. ASVspoof 2021 Workshop, 2021, pp. 61–67.
  * [10] X. Chen, Y. Zhang, G. Zhu and Z. Duan,  ``UR Channel-Robust Synthetic Speech Detection System for ASVspoof 2021,''  in Proc. ASVspoof 2021 Workshop, 2021, pp. 75–82.
  * [11] H. Tak, J. Patino, M. Todisco, A. Nautsch, N. Evans and A. Larcher,  ``End-to-end anti-spoofing with RawNet2,''  in Proc. ICASSP, 2021, pp. 6369–6373.
  * [12] W. Ge, J. Patino, M. Todisco and N. Evans,  ``Raw differentiable architecture search for speech deepfake and spoofing detection,''  in Proc. ASVspoof 2021 Workshop, 2021, pp. 22–28.
  * [13] G. Hua, A. B.-j. Teoh and H. Zhang,  ``Towards end-to-end synthetic speech detection,''  IEEE Signal Processing Letters, vol. 28, pp. 1265–1269, 2021.
  * [14] H. Tak, J.-w. Jung, J. Patino, M. Kamble, M. Todisco and N. Evans,  ``End-to-end spectro-temporal graph attention networks for speaker verification anti-spoofing and speech deepfake detection,''  in Proc. ASVspoof 2021 Workshop, 2021, pp. 1–8.
  * [15] J.-w. Jung, H.-S. Heo, H. Tak, H.-j. Shim, J. S. Chung, B.-J. Lee et al.,  ``AASIST: Audio anti-spoofing using integrated spectro-temporal graph attention networks,''  in Proc. ICASSP (to appear), 2022.
  * [16] A. B. Arrieta, N. Díaz-Rodríguez, J. Del Ser, A. Bennetot, S. Tabik, A. Barbado et al.,  ``Explainable artificial intelligence (XAI): Concepts, taxonomies, opportunities and challenges toward responsible AI,''  Information Fusion, vol. 58, pp. 82–115, 2020.
  * [17] A. Adadi and M. Berrada,  ``Peeking inside the black-box: a survey on explainable artificial intelligence (XAI),''  IEEE access, vol. 6, pp. 52138–52160, 2018.
  * [18] B. Chettri, S. Mishra, B. L. Sturm and E. Benetos,  ``Analysing the predictions of a CNN-based replay spoofing detection system,''  in 2018 IEEE Spoken Language Technology Workshop (SLT), 2018, pp. 92–97.
  * [19] Y. Zhang, W. Wang and P. Zhang,  ``The effect of silence and dual-band fusion in anti-spoofing system,''  in Proc. Interspeech, 2021, pp. 4279–4283.
  * [20] S. Becker, M. Ackermann, S. Lapuschkin, K. Müller and W. Samek,  ``Interpreting and explaining deep neural networks for classification of audio signals,''  CoRR, vol. abs/1807.03418, 2018.
  * [21] H. Muckenhirn, A. Abrol, M. Magimai-Doss and S. Marcel,  ``Understanding and visualizing raw waveform-based CNNs,''  in Proc. Interspeech, 2019, pp. 2345–2349.
  * [22] S. Sivasankaran, E. Vincent and D. Fohr,  ``Explaining deep learning models for speech enhancement,''  in Proc. Interspeech, 2021, pp. 696–700.
  * [23] S. M. Lundberg, S.-I. Lee and D. Fohr,  ``A unified approach to interpreting model predictions,''  in Advances in Neural Information Processing Systems, 2017, pp. 4765–4774.
  * [24] W. Ge, J. Patino, M. Todisco and N. Evans,  ``Explaining deep learning models for spoofing and deepfake detection with SHapley Additive exPlanations,''  in Proc. ICASSP (to appear), 2022.
  * [25] M. Christoph,  ``Interpretable machine learning: A guide for making black box models explainable,'' 2020.
  * [26] H. Dinkel, N. Chen, Y. Qian and K. Yu,  ``End-to-end spoofing detection with raw waveform CLDNNS,''  in Proc. ICASSP, 2017, pp. 4860–4864.
  * [27] Y. Ma, Z. Ren and S. Xu,  ``RW-Resnet: A novel speech anti-spoofing model using raw waveform,''  in Proc. Interspeech, 2021, pp. 4144–4148.
  * [28] Z. Teng, Q. Fu, J. White, M. E. Powell and D. C. Schmidt,  ``Complementing handcrafted features with raw waveform using a light-weight auxiliary model,''  in arXiv preprint arXiv:2109.02773, 2021.
  * [29] M. Todisco, X. Wang, V. Vestman, M. Sahidullah, H. Delgado, A. Nautsch et al.,  ``ASVspoof 2019: Future horizons in spoofed and fake audio detection,''  in Proc. Interspeech, 2019, pp. 1008–1012.
  * [30] D. Kingma and J. Ba,  ``Adam: A method for stochastic optimization,''  in Proc. International Conference on Learning Representations (ICLR), 2017, pp. 1–15.
  * [31] X. Wang and J. Yamagishi,  ``A comparative study on recent neural spoofing countermeasures for synthetic speech detection,''  in Proc. Interspeech, 2021, pp. 4259–4263.
  * [32] X. Wang, J. Yamagishi, M. Todisco, H. Delgado, A. Nautsch, N. Evans et al.,  ``ASVspoof 2019: A large-scale public database of synthesized, converted and replayed speech,''  Computer Speech & Language, vol. 64, 2020.
  * [33] A. Oord, S. Dieleman, H. Zen, K. Simonyan, O. Vinyals, A. Graves et al.,  ``Wavenet: A generative model for raw audio,''  in arXiv preprint arXiv:1609.03499, 2016.
  * [34] N. Müller, F. Dieckmann, P. Czempin, R. Canals, K. Böttinger and J. Williams,  ``Speech is silver, silence is golden: What do ASVspoof-trained models really learn?,''  in Proc. ASVspoof 2021 Workshop, 2021, pp. 55–60.
  * [35] E. Morise, F. Yokomori and K. Ozawa,  ``WORLD: A vocoder-based high-quality speech synthesis system for real-time applications,''  IEICE Trans. on Information and Systems, vol. 99, no. 7, pp. 1877–1884, 2016.
  * [36] M. Schröder, M. Charfuelan, S. Pammi and I. Steiner,  ``Open source voice creation toolkit for the MARY TTS Platform,''  in Proc. Interspeech, 2011, pp. 3253–3256.
  * [37] Z. Wu, T. Kinnunen, N. Evans, J. Yamagishi, C. Hanilçi, M. Sahidullah et al.,  ``ASVspoof 2015: the first automatic speaker verification spoofing and countermeasures challenge,''  in Proc. Interspeech, 2015, pp. 2037–2041.
  * [38] D. Matrouf, J. Bonastre and C. Fredouille,  ``Effect of speech transformation on impostor acceptance,''  in Proc. ICASSP, 2006, pp. 933–936.