# A Quasi-periodic Oscillation in the $\gamma$-ray Emission from the Non-
blazar Active Galactic Nucleus PKS 0521-36

Haiyun Zhang  Department of Astronomy, Key Laboratory of Astroparticle Physics
of Yunnan Province, Yunnan University,
Kunming 650091, China  Dahai Yan  Key Laboratory for the Structure and
Evolution of Celestial Objects, Yunnan Observatory, Chinese Academy of
Sciences,
Kunming 650011, China  Pengfei Zhang  Department of Astronomy, Key Laboratory
of Astroparticle Physics of Yunnan Province, Yunnan University,
Kunming 650091, China  Shenbang Yang  Department of Astronomy, Key Laboratory
of Astroparticle Physics of Yunnan Province, Yunnan University,
Kunming 650091, China  Li Zhang  Department of Astronomy, Key Laboratory of
Astroparticle Physics of Yunnan Province, Yunnan University,
Kunming 650091, China

###### Abstract

Long-term $\gamma$-ray variability of a non-blazar Active Galactic Nucleus
(AGN) PKS 0521-36 is investigated by using Fermi-LAT pass 8 data covering from
2008 August to 2021 March. The results show that the histogram of the
$\gamma$-ray fluxes follows a log-normal distribution. Interestingly, in the
analysis of $\sim$5.8-year (from MJD 56317 to 58447) LAT data between two
outbursts (occurring during 2012 October and 2019 May respectively), a quasi-
periodic oscillation (QPO) with a period of $\sim 1.1$ years ($\sim 5\sigma$
of significance) is found in the Lomb-Scargle Periodogram (LSP), the Weighted
Wavelet Z-transform (WWZ) and the REDFIT results. This quasi-periodic signal
also appears in the results of Gaussian process modeling the light curve.
Therefore, the robustness of the QPO is examined by four different methods.
This is the first $\gamma$-ray QPO found in a mildly beamed jet. Our results
imply that the $\gamma$-ray outbursts play an important role in the formation
of the $\gamma$-ray QPO.

Active galactic nuclei (16), Gamma-rays (637), Time series analysis (1916),
Period search (1955), Jets (870)

##  1 Introduction

It is generally accepted that Active Galactic Nuclei (AGNs) existing in the
center of a small fraction of galaxies are powerful emitters for the entire
galaxies (e.g., Padovani et al., 2017). Usually, an AGN is powered by a
central super massive black hole (SMBH) and emits variable emission covering
from radio to high-energy gamma-ray bands (e.g., Barkov et al., 2010; Madejski
& Sikora, 2016). Non-thermal GeV-TeV $\gamma$-rays from AGNs are produced in
their powerful relativistic jets.

The largest class of AGNs detected by Fermi Large Area Telescope (Fermi-LAT)
is of blazar type (Ajello et al., 2020). A typical blazar has relativistic jet
oriented at very small angle ($\theta$ $\textless 5^{\circ}$) to the line of
our sight (Ackermann et al., 2015), and therefore its jet emission is strongly
Doppler-boosted (Urry & Padovani, 1995). Rapid flares on timescale of a few
minutes in blazars have been observed at GeV $\gamma$-ray energies (e.g.,
Ackermann et al., 2016; Meyer et al., 2019) and TeV energies (e.g., Aleksić et
al., 2011). This kind of rapid flare suggests an extremely compact
$\gamma$-ray emission region which likely locates in the broad-line region
(BLR) of flat spectrum radio quasars (FSRQs). However, no signs of
$\gamma$-ray absorption by BLR photons in $\gamma$-ray spectrum suggests that
$\gamma$-ray emission region should be outside of the BLR (e.g., Costamante et
al., 2018). To reconcile the observations of fast variability and simultaneous
energy spectrum, new models are proposed (see Böttcher, 2019; Blandford et
al., 2019, for recent reviews). Moreover, timing analysis, such as the
exploration of periodicity, would provide important information about the
nature of these sources. Quasi-periodic variabilities of blazars in
$\gamma$-ray energies have been investigated, and possible quasi-periodic
variabilities in more than 20 blazars have been reported based on analyses of
Fermi-LAT data (e.g., Sandrinelli et al., 2014, 2016, 2017; Ackermann et al.,
2015; Prokhorov & Moraghan, 2017; Zhang et al., 2017a, c, b; Zhou et al.,
2018; Bhatta, 2019; Peñil et al., 2020; Zhang et al., 2020). However, the
cause for the $\gamma$-ray quasi-periodic variabilities remains controversial.
Yan et al. (2018) investigated the relation between $\gamma$-ray variability
amplitude and spectral index for PG 1553+113, and suggested that the
$\gamma$-ray quasi-periodic variability of PG 1553+113 (Ackermann et al.,
2015) could be caused by a periodic modulation in Doppler factor. Several
possible explanations have been proposed for the $\gamma$-ray quasi-periodic
variabilities in blazars, for instance, presence of a binary system of SMBHs
(Valtonen et al., 2016; Tavani et al., 2018), modulation in the accretion flow
feeding the jet (Gracia et al., 2003), and a periodic lighthouse effect
(Camenzind & Krockenberger, 1992).

Seventy non-blazar AGNs are included in the Fourth Catalog of AGNs detected by
Fermi-LAT (4LAC; Ajello et al., 2020). Among them, there are 41 radio galaxies
(RGs), 9 narrow-line seyfert 1 galaxies (NLSy1s), 2 steep-spectrum radio
quasars (SSRQs), 5 compact steep spectrum radio sources (CSSs), 1 Seyfert
galaxy, and 11 other AGNs. 6 NLSy1s, 7 RGs, 3 CSSs, 1 SSRQ, and 3 other AGNs
are found to be variable (Ajello et al., 2020). The jet properties of NLSy1
are similar to those of FSRQs, but with lower powers (e.g., D’Ammando et al.,
2013, 2015a, 2016; Zhang et al., 2013; Paliya et al., 2019). $\gamma$-ray
flares on a daily timescale are found in several NLSy1s (D’Ammando, 2019). RGs
are misaligned AGNs, and therefore the jet Doppler boosting effect is not
strong. Intraday variability is found in the Fermi-LAT data of the RG NGC 1275
(Sahakyan et al., 2018; Tanada et al., 2018). All of the non-blazar AGNs are
not bright enough to have good sampling in their long-term $\gamma$-ray light
curves, except for NGC 1275 and PKS 0521-36. The $\gamma$-ray variability and
spectral evolution of NGC 1275 have been studied in detail (e.g., Baghmanyan
et al., 2017; Sahakyan et al., 2018; Tanada et al., 2018). Here, we focus on
the $\gamma$-ray properties of PKS 0521-36. A combination of $\gamma$-ray
properties of blazars and non-blazar AGNs can provide insightful information
on AGN jet physics.

PKS 0521-36 is a peculiar non-blazar AGN with broad emission lines in the
optical and ultraviolet bands and the steep radio spectrum (Scarpa et al.,
1995; Keel, 1986). D’Ammando et al. (2015b) raised a possibility that PKS
0521-36 is an intermediate source between broad-line radio galaxy and SSRQ. It
was previously classified as an N galaxy, then a BL Lac object (D’Ammando et
al., 2015b). The latest LAT source catalog (4FGL; Abdollahi et al. 2020)
suggests that the existing data cannot determine a clear classification of
this source. Very Long Baseline Array (VLBA) observations reveal that the
radio structure of its jet is similar to misaligned AGN (D’Ammando et al.,
2015b). Recently, Angioni et al. (2019) confirmed that the jet of PKS 0521-36
is not highly beamed.

The jet structure of PKS 0521-36 has been studied widely at frequencies from
optical to X-rays (e.g., Scarpa et al., 1999; Falomo, 1994; Keel, 1986).
Because of its prominent radio, optical, and X-ray jet (Danziger et al., 1979;
Macchetto et al., 1991; Birkinshaw et al., 2002; Liuzzo et al., 2017), PKS
0521-36 is considered to be one of the most remarkable extragalactic objects.
At $\gamma$-ray energies, a strong flare on variability timescale of 12-hour
was detected by Fermi-LAT from PKS 0521-36 in 2010 June (D’Ammando et al.,
2015b). Recently, a $\gamma$-ray flare on timescale of $\sim 6$-hour in 2012
October was reported by Angioni et al. (2019).

A study on the long-term variability behavior of PKS 0521-36 in $\gamma$-ray
energies is lack. The characteristics of $\gamma$-ray variability are helpful
to understand the nature of this source. In this work, we perform a full time-
domain analysis on the latest Fermi-LAT data of PKS 0521-36. The paper is
organized as follows. The basic introduction of Fermi-LAT and results of
searching for periodicity by different methods are presented in Section 2. Our
discussion of the results is given in Section 3. Finally, a summary is
presented in Section 4.

##  2 Data Analysis and Results

###  2.1 Fermi-LAT Likelihood Analysis

Fermi-LAT is a pair conversion $\gamma$-ray detector with wide field of view
($\sim 2.4\rm\ sr$), large effective area ($>8000\ \rm cm^{2}$ at $\sim$1
GeV), and broad-energy covering from $\sim$ 20 MeV to $>$ 300 GeV. Thanks to
its all-sky monitoring capabilities, one can get long-coverage and uniform
observations of $\gamma$-ray sources. A detailed description of LAT can be
found in Atwood et al. (2009).

The data used here were collected during $\sim$12.6 years of Fermi-LAT
operation from 2008 August 4 to 2021 March 29 [Modified Julian Day (MJD):
$54682-59302$] in the 0.1-500 GeV energy range. The analysis of data follows
the standard criteria for the point-source
analysis111http://fermi.gsfc.nasa.gov/ssc/data/analysis/documentation/Pass8_usage.html.
We employ the Science Tools package of version v11r05p3 available from Fermi
Science Support
Center222https://fermi.gsfc.nasa.gov/ssc/data/analysis/software (FSSC), and
use the Pass 8 LAT database in which we choose a region of interest (ROI) of
$15^{\circ}$ centered at PKS 0521-36. Two background templates, Galactic and
extragalactic diffuse emissions, are modeled with the file gll_iem_v07.fits
and iso_P8R3_SOURCE_ V2_v1.txt, respectively. The instrumental response
function “P8R3_SOURCE_V2” is adopted in data processing. In order to get
events with high probability of being photons, we only keep the SOURCE class
events through setting “evclass = 128, evtype = 3” in the Fermi tool gtselect.
In addition, the events with zenith angle $>90^{\circ}$ are filtered out to
avoid contaminating from Earth’s limb. Good time intervals (GTIs) are selected
by using the criterion of “($\rm DATAQUAL>0$)&&($\rm LATCONFIG==1$)” in the
Fermi tool gtmktime. Above settings ensure that we adopt high-quality
observations.

###  2.2 $\gamma$-ray SED and Light Curves

In the process of analyzing the entire data, an initial binned maximum
likelihood analysis (Mattox et al., 1996) is performed. The spectral
parameters and the normalization factors of all sources lying within
$5^{\circ}$ of PKS 0521-36 as well as the normalization factors of the sources
located beyond $5^{\circ}$ but within $10^{\circ}$ from our target are set to
be free, while the sources from ROI but beyond $10^{\circ}$ have their
parameters fixed to the values in Fermi-LAT fourth source catalog (Abdollahi
et al., 2020). In the iterative fittings, the best-fitting results are
obtained, and are saved to be a new model file which is the basis of
constructing light curves and spectrum energy distribution (SED). The spectral
model of LogParabola (LP) form is used here, which is

$\frac{dN}{dE}=N_{0}\left(\frac{E}{E_{\rm
b}}\right)^{-[\alpha+\beta\log(E/E_{\rm b})]}\ .$ (1)

The global fitting results in an integrated average flux of $(1.06\pm
0.02)\times 10^{-7}$ photons cm-2 s-1 and a total Test Statistics (TS) value
of 18447.2. The TS is defined as 2log($L/L_{0}$), where $L$ is the maximum
likelihood of the model with a point source at the target position, and
$L_{0}$ is the maximum likelihood without the source (Mattox et al., 1996).
The best-fitting spectral parameters are the spectral slope $\alpha=2.43\pm
0.02$ and the curvature parameter $\beta=0.05\pm 0.01$. The average SED of PKS
0521-36 in the whole time period is shown in Figure 1.

Figure 1: Average SED of PKS 0521-36 from MJD 54682 to 59302. The best-fitting
LP model is plotted as blue solid line, and the best-fitting parameters are
given in the text. Figure 2: LAT light curve of PKS 0521-36 from MJD 54697 to
59287 in the energy range of 0.1 GeV - 500 GeV with 30 days time bin. Two
vertical dashed lines mark the time range of MJD 56317$-$58447, during which
periodic modulation is analyzed (see Section 2.4).

The unbinned maximum likelihood optimization is used to produce $\gamma$-ray
light curve of PKS 0521-36. In each time bin, the new model file mentioned
above is employed, but the spectral parameters of all sources are freezed. The
time bins having $TS<$25 and the number of the model-predicted photon $N_{\rm
pred}<$10 are excluded . This ensures that the robust detections of the source
are adopted in the analysis (Kapanadze et al., 2020; Raiteri et al., 2013). A
30-day binning light curve is produced. The light curve shown in Figure 2
reveals flux variability over the whole observational period. Three high
flaring activities appear in the light curve. The first is in 2010 June and,
has been reported by D’Ammando et al. (2015b) who also reported the
contemporaneous increases of activities observed in optical, UV and X-ray
bands. The second flaring occurred in 2012 October has been studied by Angioni
et al. (2019), and the latest flaring activity occurred in 2019 May.

Figure 3: Normalized histogram of $\gamma$-ray photon fluxes fitted by the
log-normal (red solid line) and Gaussian (blue solid line) function.

###  2.3 Flux Distribution

One important property of a time series is its probability density function
(PDF). One can explore the mechanism that drives variability by analyzing the
flux distribution of a light curve (e.g., Rieger, 2019). PDF can be estimated
by fitting a histogram of observed data. The distribution of the $\gamma$-ray
photon fluxes of PKS 0521-36 is presented in a histogram with 25 bins, and is
fitted by log-normal and Gaussian functions in the form of

$L(\phi\mid\mu,\sigma)=\frac{1}{\sqrt{2\pi}\sigma\phi}\exp\left[-\frac{(\log(\phi)-\mu)^{2}}{2\sigma^{2}}\right]\;,$
(2)
$G(\phi\mid\mu,\sigma)=\frac{1}{\sqrt{2\pi}\sigma}\exp\left[-\frac{(\phi-\mu)^{2}}{2\sigma^{2}}\right]\;,$
(3)

respectively. The parameter $\sigma$ is the standard deviation and the
parameter $\mu$ is the mean of the distribution. The flux histogram and best-
fitting results are shown in Figure 3. The errors of data points are estimated
by the method given by Gehrels (1986). The best-fitting parameters as well as
the reduced $\chi^{2}$ obtained from the fittings with Gaussian and log-normal
models are listed in table 1. It is found that the log-normal model is
preferred over the Gaussian model.

Table 1: Best-fitting parameters for the log-normal and Gaussian distribution.

Time bin | Log-normal |  | Gaussian
---|---|---|---
| Mean∗ | $\sigma^{*}$ | $\chi^{2}_{\rm red}$ |  | Mean∗ | $\sigma^{*}$ | $\chi^{2}_{\rm red}$
30-d | $2.28\pm 0.04$ | $0.49\pm 0.04$ | 0.92 |  | $9.04\pm 0.75$ | $5.27\pm 0.71$ | 1.60

Note. — ∗: in unit of $\times 10^{-8}$ photons cm-2 s-1.

###  2.4 Searching for Periodicity

Figure 4: Upper panel: $\gamma$-ray light curve from MJD 56317 to 58447 with
30-day binning. Right lower panel: two dimensional contour map of the WWZ
power spectrum of the $\gamma$-ray light curve. Left lower panel: LSP (red
dashed line) and average WWZ (black solid line) powers of the light curve.
Figure 5: LSP power calculated from the light curve during MJD 56317-58447
(blue solid line), and the 4$\sigma$ (orange dashed line), 5$\sigma$ (green
dashed line) significance curves calculated based on the simulation of
$10^{5}$ light curves with the method given by Emmanoulopoulos et al. (2013).
Figure 6: Folded light curve constructed from binned likelihood analysis of
the data from MJD 56317-58447 above 100 MeV with a period of 411 days. The red
dashed line is the constant fitting to the folded light curve in one phase.
For clarity, we show two period cycles. Figure 7: Results of the periodicity
analysis by REDFIT. The black solid line is bias-corrected power spectrum.
Dashed curves starting from the bottom are the theoretical red-noise spectrum,
90$\%$, 95$\%$ and 99$\%$ significance levels, respectively.

The Lomb-Scargle periodogram (LSP; Lomb, 1976; Scargle, 1982) and the Weighted
Wavelet Z-transform (WWZ; Foster, 1996) are two of the most employed and best
known methods for detection of periodicity in time-series in astronomy. The
LSP and WWZ methods are then applied to the light curve of PKS 0521-36 to
search for periodicity. A clear periodic modulation appears from MJD 56317 to
58447 (Figure 4). A significant peak at $\sim$400 days appears in both LSP and
WWZ powers, which indicates a quasi-periodic oscillation (QPO). The period is
estimated by fitting the power peak with a Gaussian function. It is $411\pm
33$ days given by LSP and $417\pm 39$ days given by WWZ. The bootstrap
approximation method is used to calculate the false alarm probability (FAP),
which is a robust way to estimate FAP (e.g., Prokhorov & Moraghan, 2017; Peñil
et al., 2020). The FAP is 0.0002 in the LSP power, corresponding to the 99.98%
significance level for the signal.

The algorithm given by Emmanoulopoulos et al. (2013) is also used to evaluate
the significance of this QPO signal. The properties of observed variation,
power spectral density (PSD) and PDF, are well applied in this method. The PSD
model of the form $P(f)=Af^{-\alpha}$ is used, where $\alpha$ is the power-law
index. The PDF is in the form of a log-normal distribution. The 4$\sigma$ and
5$\sigma$ significance curves are obtained by simulating the 30-day binning
light curve $10^{5}$ times, and are shown in Figure 5. The significance of the
QPO signal is $\approx 5\sigma$.

Further evidence for this QPO is presented in the folded light curve
constructed by performing the likelihood analysis method with a 411-day
period. Phase zero is set at MJD 56317. The variation with the phase is clear
(see Figure 6). The fitting to the folded light curve in one phase with a
constant results in $\chi^{2}_{\rm red}$=449.26/14, suggesting a significant
variability in the folded light curve.

An additional method, REDFIT (Schulz & Mudelsee, 2002), is also employed for
the periodicity detection. This method estimates the red-noise spectrum by
fitting the data with a first-order autoregressive (AR1) process. It can
precisely evaluate the significance of the peaks in the PSD against the red-
noise background. We use the program
REDFIT3.8e333https://www.marum.de/Prof.-Dr.-michael-schulz/Michael-Schulz-
Software.html. The bias of Fourier transform for unevenly spaced data is also
removed in the program. Two conditions should be satisfied when the program is
used for data: (1) the data points of light curve should not be too clustered;
and (2) the noise background of light curve can be approximated by an AR1
process. Obviously, our data meets the first condition. Moreover, the result
of the non-parametric runs test for our data is rtest = 15 (inside the 98$\%$
acceptance region [9,22]), which indicates that the noise background spectrum
is consistent with an AR1 model (Schulz & Mudelsee, 2002). Therefore, REDFIT
can be used to detect periodicity in our data. We set $n_{50}$=2 and select a
hanning window to reduce spectral leakage (Schulz & Stattegger, 1997). The
results show that there is a peak with more than 99$\%$ significance level,
and the corresponding periodicity is 440.9$\pm 78.6$ days (Figure 7). This
period is consistent with the LSP and the WWZ results. Note that the REDFIT
method only provides a maximum of significance of 99%.

###  2.5 Gaussian Process Modeling the $\gamma$-ray Light Curve

Besides LSP, WWZ and REDFIT, an alternative approach for the analysis of
astronomical variability is the Gaussian process modeling light curves, which
treats the variability in the time rather than frequency domain (e.g., Kelly
et al., 2009, 2014; Li & Wang, 2018). The PSD can be constructed by using
modeling results. One popular Gaussian process model is the continuous time
autoregressive moving average (CARMA) model developed by Kelly et al. (2014),
which has been applied to $\gamma$-ray variability of blazars (e.g., Goyal et
al., 2018; Ryan et al., 2019; Yang et al., 2021). Here a newly developed
Gaussian process model celerite (Foreman-Mackey et al., 2017), is used to
investigate the $\gamma$-ray variability of PKS 0521-36.

In this work, a stochastically driven damped simple harmonic oscillator (SHO)
is used, and its differential equation is given by

$\left[\frac{d^{2}}{dt^{2}}+\frac{\omega_{0}}{Q}\frac{d}{dt}+\omega_{0}^{2}\right]y(t)=\epsilon(t)\;,$
(4)

where $\omega_{0},Q,\epsilon(t)$ are the frequency of the undamped oscillator,
the quality factor of the oscillator, and a stochastic driving force (white
noise here), respectively. The PSD in the case is

$S(\omega)=\sqrt{\frac{2}{\pi}}\frac{S_{0}\omega_{0}^{4}}{(\omega^{2}-\omega_{0}^{2})^{2}+\omega^{2}\omega_{0}^{2}/Q^{2}}\;,$
(5)

where $S_{0}$ is proportional to the power when $\omega=\omega_{0}$.

A specific model for a time series can be built up with a mixture of
determined number of SHO terms. The Markov Chain Monte Carlo (MCMC) algorithm
(Foreman-Mackey et al., 2013) is used to perform fitting to light curve. From
MCMC sampling, we can evaluate values and uncertainties of model parameters.
In detail, $5\times 10^{4}$ samples are generated in our analysis, and the
first $2\times 10^{4}$ samples are taken as burn-in sampling.

Figure 8: celerite modeling results for the 30-day binning $\gamma$-ray light
curve of PKS 0521-36 from MJD 56317 to 58447; in the top plot, black points
with errors are LAT data, and the bright orange line is the best fit to the
data with the celerite model; black dots, blue histogram and orange line in
the middle plot are standardized residuals of each bin, scaled standardized
residuals distribution, and the expected normal distribution (scaled),
respectively; the ACFs of residuals and squared residuals are given in the
bottom plot where the gray region is the 95$\%$ confidence limit of the white
noise.  Figure 9: Posterior probability densities of parameters obtained from
celerite modeling results.  Figure 10: PSD constructed from the celerite
modeling results for the light curve during MJD 56317-58447. Bright color
region is the 68$\%$ confidence interval. A distinct peak appears at $\sim
400$ days.

Figure 8 shows the fitting results to the 30-day binning light curve of PKS
0521-36 from MJD 56317 to 58447 with the celerite model (SHO$\times 2$). For
the modeling results, we use Kolmogorov-Smirnov (KS) test to check the
consistency between the standardized residuals and a normal distribution with
mean zero and standard deviation of unity. The p-values of the test is p=0.23,
which means this consistency is credible (Yang et al., 2021). The auto-
correlation functions (ACFs) of the residuals and squared residuals in Figure
8 are almost inside the 95$\%$ confidence limits, so that the correlation
structures have been captured by the model. These results indicate that the
$\gamma$-ray variability can be described well by the two-SHO model. The
posterior probability densities for the model parameters are shown in Figure 9
and the value of these parameters are given in Table 2. A QPO feature with
period of $\sim 400$ days appears in the PSD constructed from the modeling
results (Figure 10).

Table 2: Parameters for the celerite modeling results corresponding to Figure
9.

Time |  SHO No. |  $lnS_{0}$ |  $lnQ$ |  $T$
---|---|---|---|---
(MJD) |  |  |  |  (days)
56317-58447 |  1 |  $-0.65_{-0.53}^{+0.54}$ |  $0.80_{-0.5}^{+0.55}$ |  $66.5_{/}^{+5.18}$
|  2 |  $-4.78_{-5.64}^{+3.57}$ |  $6.60_{-3.77}^{+5.64}$ |  $405.93_{-10.32}^{+9.68}$

Note. — (1) time period, (2) number of SHO terms, (3)-(5) parameters in
Equations 5, T=$2\pi/\omega_{0}$.

There is an aperiodic component and a periodic component in the $\gamma$-ray
variability during the time range of MJD 56317-58447. The period of the QPO is
estimated as $406\pm 10$ days which is consistent with the results given by
LSP, WWZ and REDFIT. The corresponding quality factor of the oscillator is
large (Table 2), suggesting that the QPO is significant (Foreman-Mackey et
al., 2017).

###  2.6 Short Timescale Variability

Figure 11: One-day binning light curve (blue points) and 12-hour binning light
curve (red points) for the bright flare during 2019 May 9-June 8 (MJD
58612-58642). Figure 12: Best-fitting results for the 12-hour binning light
curve showed in Figure 11 with Equation 6, and the distribution of the
residuals obtained in the fitting (bottom panel) which is fitted by a constant
(red dashed line).

Intraday $\gamma$-ray variability of PKS 0521-36 has been found in the LAT
data (D’Ammando et al., 2015b; Angioni et al., 2019). A recent bright
$\gamma$-ray flare occurred in 2019 May 9-June 8 (MJD 58612-58642) has not
been investigated. In our work, intraday variability during this flare is
observed in the one-day and half-day binning light curves (Figure 11). Here we
use unbinned maximum likelihood method to construct the light curves, and
select time bins having TS$\geqslant$ 4\. This flare is not bright enough to
use the selection criteria of TS$\geqslant$ 25 for such a short time bin.

In order to determine variability timescale, the 12-hour binning light curve
is fitted by using the following function (Abdo et al., 2010)

$F(t)=F_{\rm c}+F_{0}\left(e^{\frac{t_{0}-t}{T_{\rm
r}}}+e^{\frac{t-t_{0}}{T_{\rm d}}}\right)^{-1}\;,$ (6)

where $F_{\rm c},F_{0},t_{0},T_{\rm r}$ and $T_{\rm d}$ are a constant flux,
peak flux, peak time, and rise and decay times, respectively. This function
can be used to describe one-peak flare. The light curve with multi-peaks is
considered as a whole, and is fitted by using the model comprising several
components of Equation 6. The best-fitting results are presented in Figure 12,
and the parameters are given in Table 3. The residuals are fitted by the
constant of zero successfully with $\chi^{2}_{\rm red}=0.58$, suggesting that
our fitting to the light curve is acceptable. It is found $T_{\rm r}\sim
T_{\rm d}\sim$ 7 hours. This is the third short timescale $\gamma$-ray flare
in PKS 0521-36. However, the loose data selection criteria of TS$\geqslant$4
decreases the significance of the derived timescale.

The SEDs during 2019 May 9 - June 8 (MJD 58612-58642) and 2012 October 12 -
November 11 (MJD 56212-56242) are also produced. The spectral parameters are
listed in Table 4, together with the parameters for the average SED shown in
Figure 1. It is found that during the two flares, the SEDs are the same, which
can be described by a power-law form, different from the average LP SED.

Table 3: Best-fitting parameters for the 12-hour light curve in Figure 12.

Time | Peak No. | $F_{\rm c}$ | $F_{0}$ | $t_{0}$ | $T_{\rm r}$ | $T_{\rm d}$ | $\chi^{2}_{\rm red}$
---|---|---|---|---|---|---|---
(MJD) |  |  |  | (days) | (days) | (days) |
58620.9-58629.4 | p1 | $3.06\pm 0.54$ | $20.63\pm 8.81$ | $1.80\pm 0.31$ | $0.25\pm 0.27$ | $0.29\pm 0.18$ | 0.90
| p2 |  | $20.95\pm 5.55$ | $4.5^{*}$ | $0.24\pm 0.11$ | $0.28\pm 0.12$ |

Note. — (1) time period, (2) peak number, (3)-(7) parameters in Equations 6 ,
(8) the reduced $\chi^{2}$, (3)(4) in unit of $10^{-7}$ photons cm-2 s-1. The
value with $*$ represents that it is fixed in fitting.

Table 4: Best-fitting spectral parameters in four periods. Time period | Time period | $\alpha$ | $\beta$ | TS | Flux
---|---|---|---|---|---
(MJD) | (UT) |  |  |  | ($10^{-7}$ photons cm-2s-1)
54682-59287 | 2008-08-04/2021-03-29 | $2.43\pm 0.02$ | $0.05\pm 0.01$ | 18447.2 | $1.06\pm 0.02$
56212-56242 | 2012-10-12/2012-11-11 | $2.41\pm 0.05$ | $--$ | 1363.98 | $3.70\pm 0.21$
58612-58642 | 2019-05-09/2019-06-08 | $2.39\pm 0.08$ | $--$ | 851.68 | $3.22\pm 0.24$
56317-58447 | 2013-01-25/2018-11-25 | $2.46\pm 0.03$ | $0.06\pm 0.02$ | 6188.17 | $0.91\pm 0.02$

##  3 Discussion

Our results show that the PDF obtained from the $\gamma$-ray light curve of
PKS 0521-36 from 2008 August 4 to 2021 March 29 is preferentially log-normally
distributed. The log-normal distribution of flux is found in X-ray binaries
(e.g., Uttley & McHardy, 2001), X-ray NLSy1s (e.g., Uttley et al., 2005), and
blazars (e.g., Sinha et al., 2016; Shah et al., 2018; Ait Benkhali et al.,
2020; Bhatta & Dhital, 2020). Similar as X-ray binaries, the log-normal
distribution of $\gamma$-ray fluxes may indicate an influence of the accretion
disk on the jet (Giebels & Degrange, 2009), i.e., the multiplicative behavior
in an accretion disc impacting the jet. In this scenario, the variability
timescale is dominated by the viscous timescale of the accretion disc (Rieger,
2019). However, the tentative timescale of $\sim 7$ hours $\gamma$-ray
variability found in PKS 0521-036 is unlikely achieved in this scenario
(Czerny, 2006). The cascade-related emission processes may also lead to a log-
normal distribution of fluxes (Rieger, 2019), for instance, the cascade
emission in lepto-hadronic jet model (e.g., Mannheim, 1993; Böttcher et al.,
2013; Yan & Zhang, 2015). In this case, the observed GeV emission is produced
by the cascade processes, which is feasible. However, it would require an
extremely high jet power (e.g., Petropoulou et al., 2015; Cerruti, 2020).

A simple and natural scenario is that the log-normal distribution of
variability results from the emitting particles with the log-normal
distribution (Sinha et al., 2018; Rieger, 2019). This interpretation is
associated with the fluctuations in the particle acceleration rate (Sinha et
al., 2018). The average LP SED suggests an average LP electron distribution.
Both the log-normal and LP electron distributions could be the result of
stochastic acceleration (Tramacere et al., 2011).

Possible $\gamma$-ray QPOs have been reported in blazars (e.g., Ackermann et
al., 2015; Prokhorov & Moraghan, 2017; Peñil et al., 2020; Sandrinelli et al.,
2014, 2016, 2017; Zhang et al., 2017a, c; Zhou et al., 2018), although
cautions have been raised on the significance of these QPOs (e.g., Covino et
al., 2019, 2020; Ait Benkhali et al., 2020; Yang et al., 2021). It is well
known that blazar jets are highly beamed. So far, no $\gamma$-ray QPO is
reported in non-blazar AGNs.

Here, we report a $\gamma$-ray QPO with a significance of $\sim 5\sigma$ in
non-blazar AGN PKS 0521-36. The period is $\sim$1.1 years. This QPO behavior
only appears in the temporal range from MJD 56317 to 58447 which is between
two outbursts (occurred in 2012 October and 2019 May respectively). This
significant QPO is obtained by three different periodicity detection methods
(LSP, WWZ and REDFIT). Evidence for this QPO is also found in the Gaussian
process modeling the 30-day binning $\gamma$-ray light curve. It is confirmed
that the jet of PKS 0521-36 is a mildly beamed (Leon et al., 2016; Angioni et
al., 2019). Therefore, the QPO of PKS 0521-36 is the first $\gamma$-ray QPO
found in a mildly beamed jet. The average SED during MJD 56317-58447 is
described by a LP model (see Table 4) which is different from the power-law
SEDs in the two outbursts.

The detection of a helicoidal motion in the optical jet of PKS 0521-36 was
reported (Jiménez-Andrade et al., 2017), which may be associated with helical
magnetic field or jet precession. A simultaneous VLBI core flux activity
during the $\gamma$-ray flare in 2010-2011 was reported by Angioni et al.
(2019). This suggests that the $\gamma$-ray emission region is located inside
the radio core rather in the jet. $\gamma$-ray variability on timescale of
$\sim$7 hours is observed before and after the QPO. If the $\gamma$-ray
photons produced in a single region, the helicoidal motion along the jet at pc
scales (Jiménez-Andrade et al., 2017) could not lead to the $\gamma$-ray QPO
of PKS 0521-36.

It seems that the outburst during 2012 October plays a key role in producing
this QPO. It is not clear whether the outburst during 2019 May destroys this
QPO or triggers another QPO. The data collected in future several years will
answer this question.

##  4 Summary

We have carried out a temporal analysis of the $\gamma$-ray emissions from an
interesting non-blazar AGN PKS 0521-36 with the Fermi-LAT data from 2008
August 4 to 2021 March 29. Our main results are as follows.

$(i)$ The PDF of the $\gamma$-ray fluxes is more consistent with a log-normal
distribution rather than a Gaussian distribution. It may be due to the log-
normal distribution of emitting electrons, which could be the result of
stochastic acceleration. Moreover, the total average LP SED also can be
produced by the stochastic acceleration.

$(ii)$ Variability on timescale of $\sim 7$ hours is found during the outburst
in 2019 May. It is the third short timescale $\gamma$-ray flare in PKS
0521-36. Such a variability indicates that the $\gamma$-ray emission region is
compact. However, it should be noted that the timescale may be in low
significance.

$(iii)$ A QPO with the period of $\sim$1.1 years is found at $\sim 5\sigma$
confidence level in the time interval (MJD 56317-58447) between two outbursts.
Four different methods are used to examine the robustness of the QPO. All the
methods give significant evidence for this QPO. This is the first $\gamma$-ray
QPO found in a mildly beamed jet. The SED during the QPO interval is different
from that during the two $\gamma$-ray outbursts in 2012 October and 2019 May.
It is argued that the two $\gamma$-ray outbursts play a key role in the
formation of this $\gamma$-ray QPO. The $\gamma$-ray variability during the
QPO interval can be well modeled by two SHO components, an aperiodic component
and a periodic component.

We thank the anonymous reviewer for constructive suggestions. We acknowledge
financial support from National Key R&D Program of China under grant No.
2018YFA0404204, the National Science Foundation of China 11803081, and the
joint foundation of Department of Science and Technology of Yunnan Province
and Yunnan University [2018FY001(-003)]. The work of D. H. Yan is also
supported by the CAS Youth Innovation Promotion Association and Basic research
Program of Yunnan Province (202001AW070013). Facility: Fermi(LAT) Software:
Fermitools-conda, celerite (Foreman-Mackey et al., 2017), emcee (Foreman-
Mackey et al., 2013), REDFIT (Schulz & Mudelsee, 2002), NumPy (Harris et al.,
2020), Matplotlib (Hunter, 2007), Astropy (Astropy Collaboration et al., 2013,
2018), SciPy (Virtanen et al., 2020), DELCgen (Connolly, 2016).

## References

  * Abdo et al. (2010) Abdo, A. A., Ackermann, M., Ajello, M., et al. 2010, ApJ, 722, 520, doi: 10.1088/0004-637X/722/1/520
  * Abdollahi et al. (2020) Abdollahi, S., Acero, F., Ackermann, M., et al. 2020, ApJS, 247, 33, doi: 10.3847/1538-4365/ab6bcb
  * Ackermann et al. (2015) Ackermann, M., Ajello, M., Atwood, W. B., et al. 2015, ApJ, 810, 14, doi: 10.1088/0004-637X/810/1/14
  * Ackermann et al. (2016) Ackermann, M., Anantua, R., Asano, K., et al. 2016, ApJ, 824, L20, doi: 10.3847/2041-8205/824/2/L20
  * Ait Benkhali et al. (2020) Ait Benkhali, F., Hofmann, W., Rieger, F. M., & Chakraborty, N. 2020, A&A, 634, A120, doi: 10.1051/0004-6361/201935117
  * Ajello et al. (2020) Ajello, M., Angioni, R., Axelsson, M., et al. 2020, ApJ, 892, 105, doi: 10.3847/1538-4357/ab791e
  * Aleksić et al. (2011) Aleksić, J., Antonelli, L. A., Antoranz, P., et al. 2011, ApJ, 730, L8, doi: 10.1088/2041-8205/730/1/L8
  * Angioni et al. (2019) Angioni, R., Ros, E., Kadler, M., et al. 2019, A&A, 627, A148, doi: 10.1051/0004-6361/201935697
  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33, doi: 10.1051/0004-6361/201322068
  * Astropy Collaboration et al. (2018) Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123, doi: 10.3847/1538-3881/aabc4f
  * Atwood et al. (2009) Atwood, W. B., Abdo, A. A., Ackermann, M., et al. 2009, ApJ, 697, 1071, doi: 10.1088/0004-637X/697/2/1071
  * Baghmanyan et al. (2017) Baghmanyan, V., Gasparyan, S., & Sahakyan, N. 2017, ApJ, 848, 111, doi: 10.3847/1538-4357/aa8c7b
  * Barkov et al. (2010) Barkov, M. V., Aharonian, F. A., & Bosch-Ramon, V. 2010, ApJ, 724, 1517, doi: 10.1088/0004-637X/724/2/1517
  * Bhatta (2019) Bhatta, G. 2019, MNRAS, 487, 3990, doi: 10.1093/mnras/stz1482
  * Bhatta & Dhital (2020) Bhatta, G., & Dhital, N. 2020, ApJ, 891, 120, doi: 10.3847/1538-4357/ab7455
  * Birkinshaw et al. (2002) Birkinshaw, M., Worrall, D. M., & Hardcastle, M. J. 2002, MNRAS, 335, 142, doi: 10.1046/j.1365-8711.2002.05615.x
  * Blandford et al. (2019) Blandford, R., Meier, D., & Readhead, A. 2019, ARA&A, 57, 467, doi: 10.1146/annurev-astro-081817-051948
  * Böttcher (2019) Böttcher, M. 2019, Galaxies, 7, 20, doi: 10.3390/galaxies7010020
  * Böttcher et al. (2013) Böttcher, M., Reimer, A., Sweeney, K., & Prakash, A. 2013, ApJ, 768, 54, doi: 10.1088/0004-637X/768/1/54
  * Camenzind & Krockenberger (1992) Camenzind, M., & Krockenberger, M. 1992, A&A, 255, 59
  * Cerruti (2020) Cerruti, M. 2020, Galaxies, 8, 72, doi: 10.3390/galaxies8040072
  * Connolly (2016) Connolly, S. D. 2016, DELightcurveSimulation: Light curve simulation code.  http://ascl.net/1602.012
  * Costamante et al. (2018) Costamante, L., Cutini, S., Tosti, G., Antolini, E., & Tramacere, A. 2018, MNRAS, 477, 4749, doi: 10.1093/mnras/sty887
  * Covino et al. (2020) Covino, S., Landoni, M., Sandrinelli, A., & Treves, A. 2020, ApJ, 895, 122, doi: 10.3847/1538-4357/ab8bd4
  * Covino et al. (2019) Covino, S., Sandrinelli, A., & Treves, A. 2019, MNRAS, 482, 1270, doi: 10.1093/mnras/sty2720
  * Czerny (2006) Czerny, B. 2006, in Astronomical Society of the Pacific Conference Series, Vol. 360, AGN Variability from X-Rays to Radio Waves, ed. C. M. Gaskell, I. M. McHardy, B. M. Peterson, & S. G. Sergeev, 265
  * D’Ammando (2019) D’Ammando, F. 2019, Galaxies, 7, 87, doi: 10.3390/galaxies7040087
  * D’Ammando et al. (2013) D’Ammando, F., Orienti, M., Finke, J., et al. 2013, MNRAS, 436, 191, doi: 10.1093/mnras/stt1560
  * D’Ammando et al. (2015a) —. 2015a, MNRAS, 446, 2456, doi: 10.1093/mnras/stu2251
  * D’Ammando et al. (2015b) D’Ammando, F., Orienti, M., Tavecchio, F., et al. 2015b, MNRAS, 450, 3975, doi: 10.1093/mnras/stv909
  * D’Ammando et al. (2016) D’Ammando, F., Orienti, M., Finke, J., et al. 2016, MNRAS, 463, 4469, doi: 10.1093/mnras/stw2325
  * Danziger et al. (1979) Danziger, I. J., Fosbury, R. A. E., Goss, W. M., & Ekers, R. D. 1979, MNRAS, 188, 415, doi: 10.1093/mnras/188.2.415
  * Emmanoulopoulos et al. (2013) Emmanoulopoulos, D., McHardy, I. M., & Papadakis, I. E. 2013, MNRAS, 433, 907, doi: 10.1093/mnras/stt764
  * Falomo (1994) Falomo, R. 1994, The Messenger, 77, 49
  * Foreman-Mackey et al. (2017) Foreman-Mackey, D., Agol, E., Ambikasaran, S., & Angus, R. 2017, AJ, 154, 220, doi: 10.3847/1538-3881/aa9332
  * Foreman-Mackey et al. (2013) Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306, doi: 10.1086/670067
  * Foster (1996) Foster, G. 1996, AJ, 112, 1709, doi: 10.1086/118137
  * Gehrels (1986) Gehrels, N. 1986, ApJ, 303, 336, doi: 10.1086/164079
  * Giebels & Degrange (2009) Giebels, B., & Degrange, B. 2009, A&A, 503, 797, doi: 10.1051/0004-6361/200912303
  * Goyal et al. (2018) Goyal, A., Stawarz, Ł., Zola, S., et al. 2018, ApJ, 863, 175, doi: 10.3847/1538-4357/aad2de
  * Gracia et al. (2003) Gracia, J., Peitz, J., Keller, C., & Camenzind, M. 2003, MNRAS, 344, 468, doi: 10.1046/j.1365-8711.2003.06832.x
  * Harris et al. (2020) Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357–362, doi: 10.1038/s41586-020-2649-2
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science and Engineering, 9, 90, doi: 10.1109/MCSE.2007.55
  * Jiménez-Andrade et al. (2017) Jiménez-Andrade, E. F., Chavushyan, V., León-Tavares, J., et al. 2017, MNRAS, 470, L107, doi: 10.1093/mnrasl/slx090
  * Kapanadze et al. (2020) Kapanadze, B., Gurchumelia, A., Dorner, D., et al. 2020, ApJS, 247, 27, doi: 10.3847/1538-4365/ab6322
  * Keel (1986) Keel, W. C. 1986, ApJ, 302, 296, doi: 10.1086/163991
  * Kelly et al. (2009) Kelly, B. C., Bechtold, J., & Siemiginowska, A. 2009, ApJ, 698, 895, doi: 10.1088/0004-637X/698/1/895
  * Kelly et al. (2014) Kelly, B. C., Becker, A. C., Sobolewska, M., Siemiginowska, A., & Uttley, P. 2014, ApJ, 788, 33, doi: 10.1088/0004-637X/788/1/33
  * Leon et al. (2016) Leon, S., Cortes, P. C., Guerard, M., et al. 2016, A&A, 586, A70, doi: 10.1051/0004-6361/201527146
  * Li & Wang (2018) Li, Y.-R., & Wang, J.-M. 2018, MNRAS, 476, L55, doi: 10.1093/mnrasl/sly028
  * Liuzzo et al. (2017) Liuzzo, E., Paladino, R., & Galluzzi, V. 2017, in Submm/mm/cm QUESO Workshop 2017 (QUESO2017), 16, doi: 10.5281/zenodo.1038079
  * Lomb (1976) Lomb, N. R. 1976, Ap&SS, 39, 447, doi: 10.1007/BF00648343
  * Macchetto et al. (1991) Macchetto, F., Albrecht, R., Barbieri, C., et al. 1991, ApJ, 369, L55, doi: 10.1086/185957
  * Madejski & Sikora (2016) Madejski, G. G., & Sikora, M. 2016, ARA&A, 54, 725, doi: 10.1146/annurev-astro-081913-040044
  * Mannheim (1993) Mannheim, K. 1993, A&A, 269, 67.  https://arxiv.org/abs/astro-ph/9302006
  * Mattox et al. (1996) Mattox, J. R., Bertsch, D. L., Chiang, J., et al. 1996, ApJ, 461, 396, doi: 10.1086/177068
  * Meyer et al. (2019) Meyer, M., Scargle, J. D., & Blandford, R. D. 2019, ApJ, 877, 39, doi: 10.3847/1538-4357/ab1651
  * Padovani et al. (2017) Padovani, P., Alexander, D. M., Assef, R. J., et al. 2017, A&A Rev., 25, 2, doi: 10.1007/s00159-017-0102-9
  * Paliya et al. (2019) Paliya, V. S., Parker, M. L., Jiang, J., et al. 2019, ApJ, 872, 169, doi: 10.3847/1538-4357/ab01ce
  * Peñil et al. (2020) Peñil, P., Domínguez, A., Buson, S., et al. 2020, ApJ, 896, 134, doi: 10.3847/1538-4357/ab910d
  * Petropoulou et al. (2015) Petropoulou, M., Dimitrakoudis, S., Padovani, P., Mastichiadis, A., & Resconi, E. 2015, MNRAS, 448, 2412, doi: 10.1093/mnras/stv179
  * Prokhorov & Moraghan (2017) Prokhorov, D. A., & Moraghan, A. 2017, MNRAS, 471, 3036, doi: 10.1093/mnras/stx1742
  * Raiteri et al. (2013) Raiteri, C. M., Villata, M., D’Ammando, F., et al. 2013, MNRAS, 436, 1530, doi: 10.1093/mnras/stt1672
  * Rieger (2019) Rieger, F. 2019, Galaxies, 7, 28, doi: 10.3390/galaxies7010028
  * Ryan et al. (2019) Ryan, J. L., Siemiginowska, A., Sobolewska, M. A., & Grindlay, J. 2019, ApJ, 885, 12, doi: 10.3847/1538-4357/ab426a
  * Sahakyan et al. (2018) Sahakyan, N., Baghmanyan, V., & Zargaryan, D. 2018, A&A, 614, A6, doi: 10.1051/0004-6361/201732304
  * Sandrinelli et al. (2016) Sandrinelli, A., Covino, S., Dotti, M., & Treves, A. 2016, AJ, 151, 54, doi: 10.3847/0004-6256/151/3/54
  * Sandrinelli et al. (2014) Sandrinelli, A., Covino, S., & Treves, A. 2014, ApJ, 793, L1, doi: 10.1088/2041-8205/793/1/L1
  * Sandrinelli et al. (2017) Sandrinelli, A., Covino, S., Treves, A., et al. 2017, A&A, 600, A132, doi: 10.1051/0004-6361/201630288
  * Scargle (1982) Scargle, J. D. 1982, ApJ, 263, 835, doi: 10.1086/160554
  * Scarpa et al. (1995) Scarpa, R., Falomo, R., & Pian, E. 1995, A&A, 303, 730
  * Scarpa et al. (1999) Scarpa, R., Urry, C. M., Falomo, R., & Treves, A. 1999, ApJ, 526, 643, doi: 10.1086/308041
  * Schulz & Mudelsee (2002) Schulz, M., & Mudelsee, M. 2002, Computers and Geosciences, 28, 421, doi: 10.1016/S0098-3004(01)00044-9
  * Schulz & Stattegger (1997) Schulz, M., & Stattegger, K. 1997, Computers and Geosciences, 23, 929, doi: 10.1016/S0098-3004(97)00087-3
  * Shah et al. (2018) Shah, Z., Mankuzhiyil, N., Sinha, A., et al. 2018, Research in Astronomy and Astrophysics, 18, 141, doi: 10.1088/1674-4527/18/11/141
  * Sinha et al. (2018) Sinha, A., Khatoon, R., Misra, R., et al. 2018, MNRAS, 480, L116, doi: 10.1093/mnrasl/sly136
  * Sinha et al. (2016) Sinha, A., Shukla, A., Saha, L., et al. 2016, A&A, 591, A83, doi: 10.1051/0004-6361/201628152
  * Tanada et al. (2018) Tanada, K., Kataoka, J., Arimoto, M., et al. 2018, ApJ, 860, 74, doi: 10.3847/1538-4357/aac26b
  * Tavani et al. (2018) Tavani, M., Cavaliere, A., Munar-Adrover, P., & Argan, A. 2018, ApJ, 854, 11, doi: 10.3847/1538-4357/aaa3f4
  * Tramacere et al. (2011) Tramacere, A., Massaro, E., & Taylor, A. M. 2011, ApJ, 739, 66, doi: 10.1088/0004-637X/739/2/66
  * Urry & Padovani (1995) Urry, C. M., & Padovani, P. 1995, PASP, 107, 803, doi: 10.1086/133630
  * Uttley & McHardy (2001) Uttley, P., & McHardy, I. M. 2001, MNRAS, 323, L26, doi: 10.1046/j.1365-8711.2001.04496.x
  * Uttley et al. (2005) Uttley, P., McHardy, I. M., & Vaughan, S. 2005, MNRAS, 359, 345, doi: 10.1111/j.1365-2966.2005.08886.x
  * Valtonen et al. (2016) Valtonen, M. J., Zola, S., Ciprini, S., et al. 2016, ApJ, 819, L37, doi: 10.3847/2041-8205/819/2/L37
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261, doi: 10.1038/s41592-019-0686-2
  * Yan & Zhang (2015) Yan, D., & Zhang, L. 2015, MNRAS, 447, 2810, doi: 10.1093/mnras/stu2551
  * Yan et al. (2018) Yan, D., Zhou, J., Zhang, P., Zhu, Q., & Wang, J. 2018, ApJ, 867, 53, doi: 10.3847/1538-4357/aae48a
  * Yang et al. (2021) Yang, S., Yan, D., Zhang, P., Dai, B., & Zhang, L. 2021, ApJ, 907, 105, doi: 10.3847/1538-4357/abcbff
  * Zhang et al. (2013) Zhang, J., Liang, E.-W., Sun, X.-N., et al. 2013, ApJ, 774, L5, doi: 10.1088/2041-8205/774/1/L5
  * Zhang et al. (2017a) Zhang, P.-f., Yan, D.-h., Liao, N.-h., & Wang, J.-c. 2017a, ApJ, 835, 260, doi: 10.3847/1538-4357/835/2/260
  * Zhang et al. (2017b) Zhang, P.-f., Yan, D.-h., Liao, N.-h., et al. 2017b, ApJ, 842, 10, doi: 10.3847/1538-4357/aa7465
  * Zhang et al. (2017c) Zhang, P.-F., Yan, D.-H., Zhou, J.-N., et al. 2017c, ApJ, 845, 82, doi: 10.3847/1538-4357/aa7ecd
  * Zhang et al. (2020) Zhang, P.-f., Yan, D.-h., Zhou, J.-n., Wang, J.-c., & Zhang, L. 2020, ApJ, 891, 163, doi: 10.3847/1538-4357/ab71fe
  * Zhou et al. (2018) Zhou, J., Wang, Z., Chen, L., et al. 2018, Nature Communications, 9, 4599, doi: 10.1038/s41467-018-07103-2


# Fairness for Text Classification Tasks with
Identity Information Data Augmentation Methods

Mohit Wadhwa, Mohan Bhambhani, Ashvini Jindal, Uma Sawant, Ramanujam Madhavan
LinkedIn Corporation  (mwadhwa, mbhambhani, ajindal, usawant,
<EMAIL_ADDRESS>
###### Abstract.

Counterfactual fairness methods address the question: How would the prediction
change if the sensitive identity attributes referenced in the text instance
were different? These methods are entirely based on generating counterfactuals
for the given training and test set instances. Counterfactual instances are
commonly prepared by replacing sensitive identity terms, i.e., the identity
terms present in the instance are replaced with other identity terms that fall
under the same sensitive category. Therefore, the efficacy of these methods
depends heavily on the quality and comprehensiveness of identity pairs. In
this paper, we offer a two-step data augmentation process where (1) the former
stage consists of a novel method for preparing a comprehensive list of
identity pairs with word embeddings, and (2) the latter consists of leveraging
prepared identity pairs list to enhance the training instances by applying
three simple operations (namely identity pair replacement, identity term
blindness, and identity pair swap). We empirically show that the two-stage
augmentation process leads to diverse identity pairs and an enhanced training
set, with an improved counterfactual token-based fairness metric score on two
well-known text classification tasks.

fairness, counterfactuals, word embeddings

##  1\. Introduction

Machine learning models learn the statistical patterns from historical
instances to predict labels for future instances. These models are intended to
learn bias; for example, a toxic text classifier is expected to understand the
presence of toxic terms and assign a higher score to those terms compared to
other non-toxic terms. However, models are not intended to discriminate
between two instances based on sensitive identity terms, such as race or
gender indication terms (Bolukbasi et al., 2016; Blodgett and O’Connor, 2017;
Webster et al., 2018). Dixon et al. (2018) showed that due to the
disproportionate distribution of terms in the training data, text
classification models could unintentionally learn bias (defined as unintended
bias) specific to the identity terms.

At the instance level, a related bias issue arises when a classifier assigns
different scores to two nearly identical instances with different identity
information. Garg et al. (2019) study counterfactual fairness by considering
the question: How would the prediction change if the sensitive identity
attributes referenced in the text instance were different? To address
counterfactual fairness, it compares three approaches to address instance-
level bias and introduces a counterfactual token fairness metric. Both
measurement and mitigation methods discussed depend heavily on the
counterfactual instances generation process. Counterfactual instances are
generated by replacing identity terms present in the instance with other
identity terms that fall in the same category. For example, replacing the term
women in the sentence ”women should rule the world” with the term men to
create a counterfactual sentence ”men should rule the world”. Counterfactual
fairness methods, thus, demand information about identity pairs (such as
women-men, gay-straight), and could extensively benefit from diverse identity
pairs information and counterfactual generation mechanisms. However, Garg et
al. (2019) limits the process to a manually prepared list of 50 identity
pairs.

Figure 1. Overview of the IIDA process, (a) augments identity pairs for a
known identity pair, (b) leverages augmented identity pairs for perturbing an
example instance with blindness and replacement operations

In this paper, we present Identity Information Data Augmentation methods, IIDA
for brevity, to bolster the text classification models against unintended bias
by automatically generating diverse counterfactual instances. Given a subset
of known identity pairs (Borkan et al., 2019; Garg et al., 2019), we first
apply a simple but powerful method to augment the identity pairs with
additional multifarious identity pairs that occur in the same context. We then
use the expanded comprehensive identity pairs list to measure the
counterfactual token fairness metric. Prepared identity pairs along with a set
of instance perturbations operations like identity blindness, replacement, and
swap, are used to produce counterfactual instances that are leveraged for
training the model either in the form of augmentation or as a logit-pairing
training scheme (Garg et al., 2019; Kannan et al., 2018). In summary, the main
contributions of our paper are:

  1. (1)

Identity Pair Augmentation (Section 3.1): We propose a novel word-embeddings
based method to generate diverse identity pairs, as shown in Figure 1a.
Diverse identity pairs are leveraged for preparing counterfactual training
instances and for measuring counterfactual token fairness (Garg et al., 2019).

  2. (2)

Identity Information Perturbations (Section 3.2): To enhance the training set
with counterfactuals information, we leverage augmented identity pairs to
perform three types of text perturbation operations on training instances, as
shown in Figure 1b.

  3. (3)

Effectiveness (Section 4): We methodically evaluate IIDA on two well-
recognized unintended bias inducing text classification tasks, demonstrating
that IIDA methods help mitigate the unintended bias while sustaining the model
performance on both the tasks.

##  2\. Related work

Garg et al. (2019) propose a fairness metric, Counterfactual Token fairness
(CTF gap), for measuring model’s bias where bias is computed by taking the
absolute difference between model’s prediction score on instance and its
counterfactual, averaged over all records. Blindness, Counterfactual
Augmentation, and Counterfactual Logit Pairing methods are also discussed by
(Garg et al., 2019) for optimizing counterfactual fairness while training text
classification models. Qian et al. (2019) edit the loss function to equalize
the probabilities of male and female words to alleviate bias in word-level
language models. Wei and Zou (2019) presents data perturbation techniques for
boosting performance of text classification models. Synonym replacement,
random insertion/swap/deletion, perturbation operations are performed by Wei
and Zou (2019). Our work leverages the perturbation operations thought to
bolster fairness metric.

##  3\. IIDA

We propose a two-step data augmentation process where we first augment known
identity pairs with diverse identity pairs to make it comprehensive. We then
apply perturbation operations to enhance the training set and for the
counterfactual token fairness metric. Table 1 presents a sampled set of known
identity pairs and other possible diverse pairs.

Known Pair |  Other Possible Pairs
---|---
man-woman |  he-she, boy-girl, guy-wife, boy-person, him-girl, father-mother
blacks-whites |  black-american, latino-mexican, hispanic-african_american, mexican-african
trump-biden |  trump-obama, palin-mitt, barack-clinton, palin-hillary, trump-hillary
Table 1. Identity Pairs

###  3.1. Identity Pair Augmentation

We present the identity pair augmentation procedure in Algorithm 1. Given a
set of known identity pairs $I$, we leverage the word embedding model (Mikolov
et al., 2013; Pennington et al., 2014; Bojanowski et al., 2017), to augment
identity pairs. For each known identity pair $(i_{1},i_{2})$, we use embedding
model $M$ to extract top $top_{k1}$ most similar words for identity terms
$i_{1}$ and $i_{2}$, and form a cartesian product (line 3) of the two
extracted words set denoted by $J$. $J$ captures the set of word pairs that
possibly occur in the same context with respect to the input pair, and to
filter out the pairs that don’t occur in the same context (i.e., to remove
noisy pairs), we perform the steps presented in lines 5-9. For each pair
$(j_{1},j_{2})$, we extract top $top_{k2}$ most similar words for identity
terms $j_{1}$ and $j_{2}$ and compute a similarity score (with sets similarity
methods like Jaccard-index) between the two sets. We consider $(j_{1},j_{2})$
a valid identity pair only if the similarity score computed before is above a
selected threshold value $\tau$. Augmented identity terms are simply a set of
unique terms that are part of augmented identity pairs. We present the
parameter values in the experiment section below and show that the proposed
procedure works well in practice.

We also investigated the word embedding model with analogy operation where
given input $term1$ and $term2$, we iterate over $top_{K}$ most similar
terms/words for $term1$, and we represent each similar term by
$term1^{\prime}$. We then find analogy $term2^{\prime}$ with the following
setting:

$w(term2^{\prime})=w(term2)+w(term1^{\prime})-w(term1)$

We considered $(term1^{\prime},term2^{\prime})$ as an augmented identity pair.
We found that even a single generic term in $top_{K}$ terms could lead to
invalid out-of-context identity pairs. For example - say we consider
$(man,woman)$ as a known identity pair, then if term $`person`$ occurs in
$top_{10}$ similar words to term $`man`$, finding analogy would lead to
erroneous terms like $`individual`,`each`,`whether`,`every`$. Algorithm 1 we
propose avoids this by doing a similarity check over $top_{k2}$ similar terms
and using a similarity threshold to avoid invalid pairs. We observed that in
practice the identity pairs generated by Algorithm 1 are diverse and valid in
most cases. We share a sampled list of the identity pairs generated by
Algorithm 1 in Appendix B section.

Algorithm 1 Identity Pair Augmentation

Notations:

$\zeta$ \- function to get top-k most similar words

$\psi$ \- function to get similarity score

Input:

$I$ \- known identity pairs

$M$ \- word embedding model

$\tau$ \- similarity score threshold

$top_{k1}$ \- top-k value for augmentation

$top_{k2}$ \- top-k value for similarity check

Output: $O$ \- augmented identity pairs

Process:

1: $O$ := $set(I)$

2: $\forall(i_{1},i_{2})\in I$

3: $J$ := $\zeta(M,i_{1},top_{k1})\times\zeta(M,i_{2},top_{k1})$

4: $\forall(j_{1},j_{2})\in J$

5: $s_{1}$ := $\zeta(M,j_{1},top_{k2})$

6: $s_{2}$ := $\zeta(M,j_{2},top_{k2})$

7: if $\psi(s_{1},s_{2})\geq\tau$ do

8: $O.add((j_{1},j_{2}))$

9: endif

10: return $O$

###  3.2. Identity Information Perturbations

We here present the details to leverage augmented identity terms and pairs for
generating counterfactual instances. We additionally perform two other
operations, blindness, and swap, for instance perturbation. On each instance
in the training set, we perform the following operations:

  1. A

Identity Pair Replacement (IPR): Randomly sample $n_{r}$ identity pairs from
the augmented identity pairs and replace all the first term occurrences with
the second term of the pair.

  2. B

Identity Term Blindness (ITB): Randomly sample $n_{b}$ identity terms from the
augmented identity terms and replace all the occurrences of the sampled terms
with the blindness placeholder. Each term could have a dedicated placeholder;
for example, a gender representing term could have a placeholder
‘GENDER_TOKEN’, and similarly, ‘NAME_TOKEN’ for a public figure name.

  3. C

Identity Pair Swap (IPS): Randomly sample $n_{s}$ identity pairs from the
augmented identity pairs and swap all the sampled pairs’ occurrences.

We selected these operations based on the research by Wei and Zou (2019), and
heuristics that counterfactual and swap instances help classifier learn and
optimize diverse data while blindness operations help hide sensitive
information specific to categories with several sub-categories (e.g., public
figure names). IPR operation on an input instance $x$, leads to the generation
of diverse counterfactual instances denoted by $\phi_{r}(x)$. Similarly, ITB
operation leads to the generation of blindness instances represented by
$\phi_{b}(x)$, and the collection of instances generated by IPS operation are
characterized by $\phi_{s}(x)$. IPR and ITB operations are also considered by
(Garg et al., 2019) but with diverse identity pairs it helps to mitigate bias,
as we demonstrate in Section 4.

We further leverage $\phi_{r}$ and $\phi_{b}$ to update the logit pairing
training scheme proposed by (Garg et al., 2019) to the following objective:

(1) $\sum_{x\in X}C(f(x),y)\\\ +\lambda_{1}\sum_{x\in
X}\mathop{\mathbb{E}}_{x_{r}\sim Unif[\phi_{r}(x)]}{\mid g(x)-g(x_{r})\mid}\\\
+\lambda_{2}\sum_{x\in X}\mathop{\mathbb{E}}_{x_{b}\sim
Unif[\phi_{b}(x)]}{\mid g(x)-g(x_{b})\mid}$

where $X$, $C$, $f$, $y$ denote the training set, cost function, classifier,
and true prediction value respectively, $g$ produces a logit, and
$\lambda_{1}$, $\lambda_{2}$ are penalty hyperparameters. Counterfactual and
blinded instances individually have shown to mitigate bias (Garg et al.,
2019), and we experiment the logit pairing training scheme with a combination
of both settings. Both settings act complementary, i.e., sensitive categories
with limited list of identity pairs (for example, race sensitive category)
could be handled by counterfactual instances, and categories with exhaustive
list of identity pairs (for example, public figure names category) could be
handled by blinded instances. Swap instances $\phi_{s}$ could also be
considered with a logit pairing scheme, but we keep it out of paper scope.

##  4\. Experiments

###  4.1. Benchmark Datasets

We perform experiments on two public text classification benchmark datasets:
(1) Toxic Comment Classification Dataset 111https://www.kaggle.com/c/jigsaw-
toxic-comment-classification-challenge, (2) Hate Speech Dataset
222https://github.com/Vicomtech/hate-speech-dataset (de Gibert et al., 2018).
We combine six classes (toxic, severe_toxic, obscene, threat, insult,
identity_hate) in the toxic comment classification dataset under positive
label class ‘Toxic‘ and work with a binary classification setting. Dataset
statistics are present in Appendix A.

###  4.2. Experimental Setting

For identity pairs augmentation, we use GloVe 42B-tokens
333https://nlp.stanford.edu/projects/glove/ embedding model (Pennington et
al., 2014) with Jaccard-index for computing similarity score between two word
sets. We set $\tau=0.25$, $top_{k1}=10$, $top_{k2}=500$, and perform
lemmatization operation on the generated identity pairs to remove similar form
redundant pairs. For Hate Speech Dataset, we keep $n_{r}=50\%$, $n_{s}=10\%$
of the augmented identity pairs, and $n_{b}=10\%$ of the augmented identity
terms, for identity replacement, swap, and blindness operations respectively.
For Toxic Comment Classification Dataset, we keep $n_{r}=10\%$, $n_{s}=10\%$,
and $n_{b}=10\%$.

We work with a set of $10$ known identity pairs (shown in Appendix Section B,
Table 10), and use the proposed Algorithm 1 to augment identity pairs (shown
in Appendix Section B, Table 11).

We conduct experiments with the Convolutional Neural Network (CNN) classifier
(Kim, 2014). Previous studies on Toxic Comment Classification and Hate Speech
Classification datasets have demonstrated good performance with the CNN
classifier (Garg et al., 2019), and therefore we use it for experimentation.
We use CNN with filters of sizes 3, 4 and 5 with each having 512 filters for
all experiments. For the Hate Speech dataset, we train for 25 epochs with
early stopping criteria. For the Toxic Comment Classification dataset, we
train for 3 epochs. Training batch size across all experiments is 128. We use
categorical cross entropy as a loss function. We use Adam with a learning rate
$3\mathrm{e}{-4}$ as an optimizer. For CTF metric calculation, we generate all
possible combinations of replacements possible with the set of augmented
identity pairs.

###  4.3. Results

We use the CTF gap metric proposed by (Garg et al., 2019) as a fairness
evaluation metric where the CTF gap, for an instance, is defined as:

(2) $CTF\,gap_{\phi_{r}(x)}=\mathop{\mathbb{E}}_{x_{r}\sim
Unif[\phi_{r}(x)]}{\mid f(x)-f(x_{r})\mid}$

CTF gap over test set is the mean of CTF gap over all individual instances
(Equation 2) with all possible valid counterfactual instances. We form test
set counterfactual instances with augmented identity pairs and IPR operation.
Proposed methods could also be tested with other perturbation analysis
methods, like Prabhakaran et al. (2019), but we limit the paper’s scope to CTF
metric.

Table 2 and 3 show performance metric and fairness metric results
respectively, for Toxic Comment Classification task. We denote Known Identity
Pairs set by ‘KIP’ and Augmented Identity Pairs set by ‘AIP’. We measure
Accuracy (at threshold = $0.5$) and CTF gap at the label class level (positive
class denoted by ‘pos’, negative class denoted by ‘neg’, and combined denoted
by ‘all’) to better understand the impact of mitigation methods. Blindness
method (CNN AIP Blindness) by design results in a CTF gap value close to 0.
Blindness method, however, misses on a lot of information that may be useful
for classification tasks, and hence, is not usually recommended for bias
mitigation. Other shortcomings of the blindness method are discussed further
in (Garg et al., 2019). The augmentation method (CNN AIP Augmentation) adds
$\phi_{r}$, $\phi_{b}$, and $\phi_{s}$ instances for each training instance to
the training set and trains the model with the augmented set. CNN AIP
Augmentation method, as observed empirically, reduces the CTF gap value. The
logit pairing method (CNN AIP LP), as expected, performs best with a
significant reduction in the CTF gap value. We present the results with three
sets of penalty values, i.e., (1) $\lambda_{1}=1$, $\lambda_{2}=0$, (2)
$\lambda_{1}=0$, $\lambda_{2}=1$, (3) $\lambda_{1}=1$, $\lambda_{2}=1$. CNN
AIP LP $\lambda_{1}=1$, $\lambda_{2}=1$ works better than other methods.
Appendix C presents the details about other values of
$\lambda_{1},\lambda_{2}$. We also show with two classifier settings (namely
CNN KIP Augmentation and CNN KIP LP) that the proposed identity pairs
augmentation method leads to improved CTF gap value. Overall, all the
mitigation methods reduce the CTF gap value without harming the classifier’s
overall Accuracy and AUC score. We observe similar results for Hate Speech
Classification task, as shown in Table 4 and 5, where CNN AIP LP
$\lambda_{1}=1$, $\lambda_{2}=1$ works best, i.e., shows low CTF gap value and
doesn’t affect model’s performance metrics. We visualize AUC vs. CTFall
comparison in Figure 2 and 3 to better understand the performance vs. fairness
metric trade-off.

Model Metric | Accuracy |
---|---|---
$Acc_{pos}$ |  $Acc_{neg}$ |  $Acc_{all}$ | AUC
CNN | 0.8271 | 0.9264 | 0.9167 | 0.9587
CNN KIP Augmentation | 0.8708 | 0.9050 | 0.9016 | 0.9584
CNN KIP LP $\lambda_{1}=1,\lambda_{2}=1$ | 0.7533 | 0.9468 | 0.9278 | 0.9534
CNN AIP Blindness | 0.8258 | 0.9217 | 0.9123 | 0.9558
CNN AIP Augmentation | 0.8520 | 0.9139 | 0.9078 | 0.9570
CNN AIP LP $\lambda_{1}=1,\lambda_{2}=0$ | 0.7804 | 0.9364 | 0.9211 | 0.9537
CNN AIP LP $\lambda_{1}=0,\lambda_{2}=1$ | 0.8117 | 0.9246 | 0.9135 | 0.9538
CNN AIP LP $\lambda_{1}=1,\lambda_{2}=1$ | 0.7327 | 0.9460 | 0.9251 | 0.9489
Table 2. Toxic Comment Classification Task : Performance Metric Results Model Metric | CTF Gap
---|---
$CTF_{pos}$ |  $CTF_{neg}$ |  $CTF_{all}$
CNN | 0.0229 | 0.0069 | 0.0084
CNN KIP Augmentation | 0.0198 | 0.0074 | 0.0086
CNN KIP LP $\lambda_{1}=1,\lambda_{2}=1$ | 0.0205 | 0.0042 | 0.0057
CNN AIP Blindness | 0.0000 | 0.0000 | 0.0000
CNN AIP Augmentation | 0.0119 | 0.0037 | 0.0045
CNN AIP LP $\lambda_{1}=1,\lambda_{2}=0$ | 0.0059 | 0.0011 | 0.0016
CNN AIP LP $\lambda_{1}=0,\lambda_{2}=1$ | 0.0064 | 0.0017 | 0.0021
CNN AIP LP $\lambda_{1}=1,\lambda_{2}=1$ | 0.0027 | 0.0007 | 0.0008
Table 3. Toxic Comment Classification Task : Fairness Metric Results Model Metric | Accuracy |
---|---|---
$Acc_{pos}$ |  $Acc_{neg}$ |  $Acc_{all}$ | AUC
CNN | 0.8368 | 0.6653 | 0.7510 | 0.8280
CNN KIP Augmentation | 0.7657 | 0.7406 | 0.7531 | 0.8159
CNN KIP LP $\lambda_{1}=1,\lambda_{2}=1$ | 0.7406 | 0.7155 | 0.7280 | 0.8156
CNN AIP Blindness | 0.7322 | 0.7490 | 0.7406 | 0.8129
CNN AIP Augmentation | 0.8075 | 0.6569 | 0.7322 | 0.8063
CNN AIP LP $\lambda_{1}=1,\lambda_{2}=0$ | 0.8159 | 0.6778 | 0.7469 | 0.8180
CNN AIP LP $\lambda_{1}=0,\lambda_{2}=1$ | 0.8075 | 0.6987 | 0.7531 | 0.8176
CNN AIP LP $\lambda_{1}=1,\lambda_{2}=1$ | 0.7874 | 0.6769 | 0.7322 | 0.8111
Table 4. Hate Speech Classification Task : Performance Metric Results Model Metric | CTF Gap
---|---
$CTF_{pos}$ |  $CTF_{neg}$ |  $CTF_{all}$
CNN | 0.0598 | 0.0776 | 0.0671
CNN KIP Augmentation | 0.0507 | 0.0668 | 0.0573
CNN KIP LP $\lambda_{1}=1,\lambda_{2}=1$ | 0.0428 | 0.0529 | 0.0469
CNN AIP Blindness | 0.0001 | 0.0000 | 0.0001
CNN AIP Augmentation | 0.0342 | 0.0264 | 0.0310
CNN AIP LP $\lambda_{1}=1,\lambda_{2}=0$ | 0.0139 | 0.0230 | 0.0177
CNN AIP LP $\lambda_{1}=0,\lambda_{2}=1$ | 0.0170 | 0.0250 | 0.0203
CNN AIP LP $\lambda_{1}=1,\lambda_{2}=1$ | 0.0110 | 0.0141 | 0.0123
Table 5. Hate Speech Classification Task : Fairness Metric Results Figure 2.
AUC vs. CTFall comparison on Toxic Comment Classification Task Figure 3. AUC
vs. CTFall comparison on Hate Speech Classification Task

##  5\. Conclusion

We demonstrate the efficacy of proposed Identity Information Data Augmentation
methods, IIDA, to help mitigate bias in text classification tasks. The two-
step augmentation setup we propose leads to diverse identity pairs and
enhanced training instances based on simple operations like identity
replacement, blindness, and swap. We further show that amalgam of
counterfactual and blindness based logit-pairing training scheme results in a
more fair classifier. We believe that IIDA’s simple methods could be
efficiently leveraged in practice. Continued work on this topic could examine
identity pairs validation, perform analysis at different subgroups level,
study the effects of out-of-vocabulary tokens, explore perturbation operations
enhancement to create and validate meaningful instances.

## References

  * (1)
  * Blodgett and O’Connor (2017) Su Lin Blodgett and Brendan O’Connor. 2017.  Racial disparity in natural language processing: A case study of social media african-american english.  _Talk at the 2017 Workshop on Fairness, Accountability, and Transparency in Machine Learning_ (2017).
  * Bojanowski et al. (2017) Piotr Bojanowski, Edouard Grave, Armand Joulin, and Tomas Mikolov. 2017.  Enriching word vectors with subword information.  _Transactions of the Association for Computational Linguistics_ 5 (2017), 135–146.
  * Bolukbasi et al. (2016) Tolga Bolukbasi, Kai-Wei Chang, James Zou, Venkatesh Saligrama, and Adam Kalai. 2016.  Man is to computer programmer as woman is to homemaker? debiasing word embeddings.  _NeurIPS_ (2016).
  * Borkan et al. (2019) Daniel Borkan, Lucas Dixon, Jeffrey Sorensen, Nithum Thain, and Lucy Vasserman. 2019.  Nuanced metrics for measuring unintended bias with real data for text classification. In _Companion Proceedings of The 2019 World Wide Web Conference_. 491–500.
  * de Gibert et al. (2018) Ona de Gibert, Naiara Perez, Aitor García-Pablos, and Montse Cuadros. 2018.  Hate speech dataset from a white supremacy forum.  _Proceedings of the 2nd Workshop on Abusive Language Online (ALW2)_ (2018).
  * Dixon et al. (2018) Lucas Dixon, John Li, Jeffrey Sorensen, Nithum Thain, and Lucy Vasserman. 2018.  Measuring and mitigating unintended bias in text classification. In _Proceedings of the 2018 AAAI/ACM Conference on AI, Ethics, and Society_. 67–73.
  * Garg et al. (2019) Sahaj Garg, Vincent Perot, Nicole Limtiaco, Ankur Taly, Ed H Chi, and Alex Beutel. 2019\.  Counterfactual fairness in text classification through robustness. In _Proceedings of the 2019 AAAI/ACM Conference on AI, Ethics, and Society_. 219–226.
  * Kannan et al. (2018) Harini Kannan, Alexey Kurakin, and Ian Goodfellow. 2018\.  Adversarial logit pairing.  _CoRR abs/1803.06373_ (2018).
  * Kim (2014) Yoon Kim. 2014.  Convolutional Neural Networks for Sentence Classification. In _Proceedings of the 2014 Conference on Empirical Methods in Natural Language Processing (EMNLP)_. Association for Computational Linguistics, Doha, Qatar, 1746–1751.  https://doi.org/10.3115/v1/D14-1181
  * Mikolov et al. (2013) Tomas Mikolov, Kai Chen, Greg Corrado, and Jeffrey Dean. 2013\.  Efficient estimation of word representations in vector space.  _International Conference on Learning Representations_ (2013).
  * Pennington et al. (2014) Jeffrey Pennington, Richard Socher, and Christopher D Manning. 2014.  Glove: Global vectors for word representation. In _Proceedings of the 2014 conference on empirical methods in natural language processing (EMNLP)_. 1532–1543.
  * Prabhakaran et al. (2019) Vinodkumar Prabhakaran, Ben Hutchinson, and Margaret Mitchell. 2019\.  Perturbation sensitivity analysis to detect unintended model biases.  _Proceedings of the Conference on Empirical Methods in Natural Language Processing (EMNLP)_ (2019).
  * Qian et al. (2019) Yusu Qian, Urwa Muaz, Ben Zhang, and Jae Won Hyun. 2019\.  Reducing gender bias in word-level language models with a gender-equalizing loss function.  _arXiv preprint arXiv:1905.12801_ (2019).
  * Webster et al. (2018) Kellie Webster, Marta Recasens, Vera Axelrod, and Jason Baldridge. 2018.  Mind the gap: A balanced corpus of gendered ambiguous pronouns.  _Transactions of the Association for Computational Linguistics_ 6 (2018), 605–617.
  * Wei and Zou (2019) Jason Wei and Kai Zou. 2019\.  Eda: Easy data augmentation techniques for boosting performance on text classification tasks.  _Proceedings of the 2019 Conference on Empirical Methods in Natural Language Processing and the 9th International Joint Conference on Natural Language Processing (EMNLP-IJCNLP)_ (2019).

##  Appendix A Experiment Details

Dataset statistics are present in Table 6 and 7.

| Toxic | Non-toxic | Total
---|---|---|---
Train | 16225 | 143346 | 159571
Test | 6243 | 57735 | 63978
Table 6. Toxic Comment Classification Dataset | Hate | No Hate | Total
---|---|---|---
Train | 957 | 957 | 1914
Test | 239 | 239 | 478
Table 7. Hate Speech Dataset

We present results on the Hate Speech Dataset for other configurations of
$\lambda_{1},\lambda_{2}$ in Table 8 and 9.

Model Metric | Accuracy |
---|---|---
$Acc_{pos}$ |  $Acc_{neg}$ |  $Acc_{all}$ | AUC
CNN AIP LP $\lambda_{1}=0.5,\lambda_{2}=0.5$ | 0.8243 | 0.6611 | 0.7427 | 0.8120
CNN AIP LP $\lambda_{1}=5,\lambda_{2}=5$ | 0.6946 | 0.5105 | 0.6025 | 0.6536
Table 8. Hate Speech Classification Task : Performance Metric Results Model Metric | CTF Gap
---|---
$CTF_{pos}$ |  $CTF_{neg}$ |  $CTF_{all}$
CNN AIP LP $\lambda_{1}=0.5,\lambda_{2}=0.5$ | 0.0151 | 0.0219 | 0.0179
CNN AIP LP $\lambda_{1}=5,\lambda_{2}=5$ | 0.0006 | 0.0007 | 0.0006
Table 9. Hate Speech Classification Task : Fairness Metric Results

##  Appendix B Generated Identity Pairs Set

Term 1 | Term 2
---|---
man | woman
democrat | republican
trump | biden
blacks | whites
asian | latino
china | america
africa | europe
christian | jewish
buddhist | atheist
homosexual | heterosexual
Table 10. Known Identity Pairs

We here present a small sample of pairs generated by Algorithm 1.

Term 1 | Term 2
---|---
female | male
religious | christianity
mexican | hispanic
she | men
japanese | asian
woman | man
she | he
muslim | christian
asian | latino
barack | mitt
gingrich | hillary
american | african-american
woman | him
liberal | conservative
american | black
girl | boy
lady | guy
gay | lesbian
jew | christian
atheist | buddhist
asia | australia
america | japan
latino | african-american
she | person
homosexual | homosexuality
lady | boy
she | him
africa | zimbabwe
australia | italy
australia | america
korea | europe
Table 11. Generated Identity Pairs