# Single Production of Vectorlike Y Quarks at the HL-LHC

V. Cetinkaya  Department of Physics, Kutahya Dumlupinar University, 43100
Kutahya, Turkey  A. Ozansoy  Department of Physics, Ankara University, 06100
Ankara, Turkey  V. Ari  Department of Physics, Ankara University, 06100
Ankara, Turkey  O. M. Ozsimsek  Department of Physics Engineering, Hacettepe
University, 06800 Ankara, Turkey  O. Cakir  Department of Physics, Ankara
University; 06100 Ankara, Turkey

###### Abstract

We study single production of exotic vectorlike $Y$ quark with electric charge
$|Q_{Y}|=4/3$ and its subsequent decay at the High Luminosity LHC (HL-LHC).
Most of the vector like quark (VLQ) decays have the electroweak $W$ bosons in
the intermediate state. Besides their direct productions singly or pairs, the
$W$-bosons are involved in decay chains as a result of the decay of a top
quark which contributes to the background. This is particularly the case since
vectorlike $Y$ quark, which is estimated to be produced with a high cross-
section, can only decay via a $W$ boson and a down type quark ($d,s,b$). We
calculate the cross sections of signal (for different couplings and mass
values) and relevant Standard Model (SM) backgrounds. After a fast simulation
of signal and background events, estimations of the sensitivity to the
parameters (mass range 1000-2500 GeV for coupling value $\kappa_{Y}=0.5$, and
mass range 500-2000 GeV for coupling values $\kappa_{Y}=0.3$ and
$\kappa_{Y}=0.15$) have been presented at the HL-LHC with center of mass
energy $\sqrt{s}=14$ TeV and integrated luminosity projections of 300 fb-1,
1000 fb-1 and 3000 fb-1.

vectorlike, Quarks, HL-LHC

###### pacs:

14.65.Jk–Other quarks, 13.85.Rm–Limits on production of particles,
12.15.Ff–Quark and lepton masses and mixing.

##  I Introduction

The results from the experiments at the Large Hadron Collider (LHC) have
confirmed the validity of the Standard Model (SM) of particle physics up to a
high energy scale and intensity. The upgrade of the LHC to the high luminosity
phase (HL-LHC) (key-1, ) at center-of-mass energy of 14 TeV and integrated
luminosity of 3000 fb-1 will extend the sensitivity and perspectives, with the
upgraded detectors and large data, to possible opportunities beyond the SM.
This upgrade will be crucial for precision measurements in the Higgs sector
and for increasing the new physics discovery potential at the energy and
luminosity frontier. The large data will be collected in two steps. During the
first upgrade phase (Run-3), experiments are expected to collect an integrated
luminosity of 300 fb-1, whereas in the next phase (Run-4), a total amount of
data corresponding to 3000 fb-1 is foreseen.

Predictions for the existence of new fermionic resonances referred to as
vectorlike quarks, which are also common in some beyond the Standard Model
(BSM) scenarios, have been expressed recently. Vectorlike quarks (key-2, ;
key-3, ) are defined as colour-triplet under $SU(3)$ and spin-1/2 fermions
whose left-handed and right-handed chiral components have the same
transformation properties under the $SU(2)\times U(1)$ gauge group.

The ATLAS (key-4, ) and CMS (key-5, ) Collaborations have published searches
for single production of vectorlike $T/Y$ quarks in decay channel $T/Y\to Wb$
and set $95\%$ confidence level (C.L.) lower limits on $T$/$Y$ quark masses.
The upper limits on the couplings are $|\sin\theta_{L}|=0.18$ for a singlet
$T$ quark, and $|\sin\theta_{R}|=0.17$ for a $(B,Y)$ doublet model, and
$|\sin\theta_{L}|=0.16$ for a $(T,B,Y)$ triplet model for a $Y$ quark mass of
$800$ GeV. Within the $(B,Y)$ doublet model, the limits on the mixing
parameter is comparable with the exclusion limits from electroweak precision
observables in the mass range [$900-1250$] GeV (key-4, ). Upper limits are
placed on the production cross section of heavy exotic quarks by the CMS
experiment (key-5, ), for $Y$ quarks with coupling of $0.5$ and $B(Y\to
bW)=100\%,$ the observed (expected) lower mass limits are given as $1400$
($1000$) GeV.

Potential of the HL-LHC in searching for the bounds on mass and coupling would
be higher. In the analysis, we investigate a resonance particle as vectorlike
$Y$ quark in the invariant mass distributions. In order to measure its charge,
leptonic decay of $W$ boson and the charge of $b$-tagged jet can be used. The
single production of vector like quarks is model dependent, the framework of
the model suggests vectorlike $Y(-4/3)$ quark to $Wb$ channel with a branching
ratio of $100\%$, however vectorlike $T(2/3)$ quark decay into the same
channel but with different branching (for example $50\%$ for singlet
configuration). This will effect the cross section times branching ratio and
one expect stronger bounds for the $Y$. Usually experimental searches for
vectorlike quarks adopt a phenomenological approach, assuming that only one
new VLQ state is present beyond the SM. The analysis has been designed based
on a simplified scenario for modelling the VLQ dynamics, assuming that only
one new VLQ is present beyond the SM for parametrising its single production.
Single production of vectorlike quarks with large width effects have already
been studied in Refs. (key-6, ; key-7, ).

In this study, we consider an effective model framework for the single
production of vectorlike $Y$ quark in the second section. In the third
section, we mention about the decay width of vectorlike $Y$ quark. Production
cross sections for signal process as well as corresponding SM backgrounds are
given in the fourth section. Modeling of the signal and background events are
performed in the fifth section. After detector simulation, event selection and
analysis results have been presented in the sixth section. Finally,
statistical significance of the signal have been given depending on the
parameter space (mass and coupling) of the model framework, and we draw a
conclusion on the search potential for vectorlike $Y$ quark at the HL-LHC.

##  II Model Framework

Depending on the model framework, vectorlike quarks are classified as $SU(2)$
singlets, doublets or triplets of flavours $B,T,X$ or $Y$, in which the first
two have the same charge as the standard model (SM) bottom and top quarks
while the vectorlike quark $X$ and $Y$ have exotic electric charge $5e/3$ and
$-4e/3$, respectively. In this framework, vectorlike quark $Y$ can exists as
($B,Y$) doublet or ($T,B,Y$) triplet. As in the description detailed in Ref.
(key-8, ), an exotic vectorlike $Y$ quark can decay into a $W$ boson and a
down sector ($d,s,b$) quark. Within the framework, vectorlike quarks are
expected to couple preferentially to third-generation quarks and can have
flavour-changing neutral-current decays in addition to the charged-current
decays characteristic of chiral quarks (key-8, ). We use an effective
Lagrangian framework for the interactions of vectorlike $Y$ quark with the SM
quarks through the $W$ boson exchange including free parameters
($\kappa_{qL}^{Y},\kappa_{qR}^{Y}$):

$\displaystyle L_{Y}$
$\displaystyle=i\bar{Y}\\!\\!\cancel{D}Y\text{\textminus}m_{Y}\bar{Y}Y$
$\displaystyle+\frac{g}{\sqrt{2}}[\bar{Y}\\!\\!\cancel{W}(\kappa_{dL}^{Y}P_{L}+\kappa_{dR}^{Y}P_{R})d+\text{h.c.}$
$\displaystyle+\frac{g}{\sqrt{2}}[\bar{Y}\\!\\!\cancel{W}(\kappa_{sL}^{Y}P_{L}+\kappa_{sR}^{Y}P_{R})s+\text{h.c.}$
$\displaystyle+\frac{g}{\sqrt{2}}[\bar{Y}\\!\\!\cancel{W}(\kappa_{bL}^{Y}P_{L}+\kappa_{bR}^{Y}P_{R})b+\text{h.c.}$
(1)

where $\cancel{D}$ and $\cancel{W}$ are the definitions $\gamma^{\mu}D_{\mu}$
and $\gamma^{\mu}W_{\mu}$ for the covariant derivative and $W$ boson field,
respectively (key-9, ). The terms in the first line denote gauge invariant
kinetic and mass term for the vectorlike $Y$ quark fields, both strong and
electroweak pieces of the covariant derivative can be included in this part.
However, the electroweak pieces are expected to yield a small effect with
respect to their strong interaction part (key-8, ). The next three lines
define the interactions between vectorlike $Y$ quark and SM down type quarks
($d,s,b$) through a $W$ boson exchange. Vectorlike quarks can also induce a
mixing between the SM and new physics (NP) sectors. The corresponding elements
of the mixing matrices can be inserted in the interactions strength. The
$\kappa_{qL}^{Y}$ and $\kappa_{qR}^{Y}$ parameters include the relevant
elements of the quark mixing matrices with the corresponding projection
operators $P_{L}$ and $P_{R}$. The mixing of VLQs simultaneously with more
than one SM family is strongly constrained by various flavor changing
processes and we focus on a mixing with the third family in the work .
Following this framework interaction terms in the Lagrangian (Eq. 1) have been
supplemented into the SM Lagrangian.

##  III Decay of VectorLike Y Quark

The decay width of vectorlike $Y$ quark has been expressed as

$\Gamma(Y\to
Wq)=\frac{\alpha_{e}(\kappa_{qL}^{2}+\kappa_{qR}^{2})}{16s_{W}^{2}}\frac{(m_{W}^{2}-m_{Y}^{2})^{2}(2m_{W}^{2}+m_{Y}^{2})}{m_{W}^{2}m_{Y}^{3}}$
(2)

where down-type SM quarks ($q=d,s,b$) are taken massless. The interaction
strength can be parametrized in terms of electromagnetic coupling constant
$g_{e}=\sqrt{4\pi\alpha_{e}}$ and $s_{W}$ (which stands for the sine of the
electroweak mixing angle). The relative importance of a decay channel, for
instance $\Gamma(Y\to Wq)$, can be expressed by the ratio $\Gamma(Y\to
Wq)/\sum_{q}\Gamma(Y\to
Wq)=(\kappa_{qL}^{2}+\kappa_{qR}^{2})/\sum_{q}(\kappa_{qL}^{2}+\kappa_{qR}^{2})$.
Taking into account the current limits on the mass and couplings of vectorlike
quarks from high energy experiments, and sensitivity to the new physics
parameters, we focus on that vectorlike $Y$ quark couples only to the bottom
quark through the charged current interacton ($\kappa_{Y}=\kappa_{bL}$ or
$\kappa_{bR}$). For large mass values of vectorlike $Y$ quark, i.e. ratio
$m_{W}/m_{Y}\ll 1$, the decay width is approximated as $\Gamma(Y\to Wb)\sim
A\kappa_{Y}^{2}m_{Y}^{3}$, where $A$ is a constant with a value of $3.28\times
10^{-7}$ GeV-2.

The decay width for vectorlike $Y$ quark depending on the mass $m_{Y}$ for
different value of the coupling parameter $\kappa_{Y}=0.5,0.2,0.1,0.05$ has
been presented in Fig. 1. As it can be seen from this figure, decay width
changes quadratically with the values of coupling parameters and changes as a
function of mass approximately ($\approx m_{Y}^{3}$) as given in Eq. 2. Given
the specified width, the couplings are implicit in the model for fixed mass
value.

Figure 1: Decay width for vectorlike $Y$ quark depending on the mass $m_{Y}$
for different value of the coupling parameter $\kappa_{Y}$.

##  IV Production Cross Section

In order to make a prediction for the signal, we calculate cross section for
on-shell vectorlike $Y$ quark production. Here, we consider $pp\to Ybj+X$
signal process for investigating interactions between vectorlike $Y$ quark and
SM bottom quark via $W$ boson specified in Eq. 1. The representative diagram
for the subprocess $qg\to Y\bar{b}q^{\prime}$ (a similar diagram for $\bar{Y}$
single production) with subsequent decays $Y\to W^{-}b\to q\bar{q}^{\prime}b$
is presented in Fig. 2. It is representative for the brevity, because there
are a lot of signal diagrams for both $Y\bar{b}j$ and $\bar{Y}bj$ productions,
corresponding to the $W^{\pm}$ exchange, light quark flavors $q$ and
$q^{\prime}$, $b$-quark and anti $b$-quark, which we take them into account in
the event generation.

Figure 2: Representative diagram for single production of vectorlike $Y$ quark
with subsequent hadronic decay channel.

The cross section for single production of vectorlike $Y$ quark through the
process $pp\to Ybj+X$ depending on different mass $m_{Y}$ and coupling values
$\kappa_{Y}$ has been presented in Fig. 3. As it can be seen from this figure,
the cross section has large values especially in the low mass region and it
increases depending on the increase in coupling parameters. In more detail, we
present the signal cross section for different mass
$m_{Y}=1000,1500,2000,2500$ GeV and coupling parameter values
$\kappa_{Y}=0.5,0.2,0.1,0.05$ in Table A.1. We present cross section numerical
values in the Appendix A to avoid similar information as in Fig. 3. The cross
section decreases according to the decreasing values of coupling parameter
$\kappa_{Y}$ and to the increasing values of the mass $m_{Y}$. While making
these cross section calculations, we take into account automatic calculation
(auto) mode of decay width $\Gamma_{Y}$. In the single production of
vectorlike $Y$ quark, which is model dependent, the cross section can be
related to a function of coupling $\kappa_{Y}$ and mass $m_{Y}$ within the
model framework.

Figure 3: Cross section for the process $pp\to Ybj+X$ depending on the mass of
vectorlike $Y$ quark for different coupling values.

##  V Modeling of Signal and Background

The signal samples are generated with the MadGraph5_aMC@NLO (key-10, ) using
VLQ_V4_4FNS_UFO model implemented in FeynRules (key-11, ; key-12, ), through
single production $pp\to Ybj$ as a high mass on-shell particle $Y$ with left-
handed or right-handed couplings and subsequent decays $Y\to Wb$ followed by
$W$ boson decaying hadronically. However, the kinematics of the final-state
jets are similar for left-handed ($\kappa_{bL}^{Y}$) and right-handed
($\kappa_{bR}^{Y}$) couplings, then the acceptances for two chiralities are
found to be the same. All hadronic channel is optimised to search for massive
vectorlike $Y$ quark which decays to a $W$ boson and a high-momentum $b$-jet
in the final state. The $W$-boson and the $b$-jet originating from the
vectorlike $Y$ quark decay are expected to be back-to-back in the transverse
plane. The signal topology includes an outgoing light quark in the process
which often produces a forward jet in the detector. The second $b$-jet in the
signal events comes from the gluon to a pair of $b$-quark ($b\bar{b}$)
splitting which may be observed in either the forward or central region.
Having typically low momentum, it is often assumed to be outside of the
detector acceptance.

Signal and background cross sections are given in Table 1. The SM background
simulation samples include top quark pair $t\bar{t}$, $W+jets$, $Z+jets$,
single top quark ($tj,$ $tb$, $tW$), associate top and $Z$ production ($tZj$),
$WW$, $ZZ$ and $WZ$ dibosons. The events for these processes are generated
with MadGraph5_aMC@NLO (key-10, ). Monte Carlo event simulations for SM
background and signal are interfaced with Pythia8 (key-13, ) for fragmentation
and showering. Signal and background samples use the NNPDF2.3 PDF set (key-14,
). Both the signal and background events are produced with the generator level
cuts: such as minimum transverse momentum of jets $p_{T}(j)=20$ GeV, maximum
pseudo-rapidity for jets $|\eta(j)|=5$, and minimum distance between jets
$\Delta R(j,j)=0.4$.

Table 1: Process, mode and generated jets in modeling the signals and different backgrounds, and calculated cross sections. While calculating these cross sections, the $W/Z$ boson and top quark are allowed to decay hadronically.  Process | Mode | Gen jets | Cross Section (pb)
---|---|---|---
$\kappa_{Y}=0.5$, $m_{Y}=1000\text{ GeV}$ | $Ybj$ | $2b+3j$ | $6.499\times 10^{-1}$
$\kappa_{Y}=0.5$, $m_{Y}=1500\text{ GeV}$ | $Ybj$ | $2b+3j$ | $1.520\times 10^{-1}$
$\kappa_{Y}=0.5$, $m_{Y}=2000\text{ GeV}$ | $Ybj$ | $2b+3j$ | $4.160\times 10^{-2}$
$\kappa_{Y}=0.5$, $m_{Y}=2500\text{ GeV}$ | $Ybj$ | $2b+3j$ | $1.236\times 10^{-2}$
$\kappa_{Y}=0.3$, $m_{Y}=500\text{ GeV}$ | $Ybj$ | $2b+3j$ | $1.422\times 10^{0}$
$\kappa_{Y}=0.3$, $m_{Y}=1000\text{ GeV}$ | $Ybj$ | $2b+3j$ | $2.367\times 10^{-1}$
$\kappa_{Y}=0.3$, $m_{Y}=1500\text{ GeV}$ | $Ybj$ | $2b+3j$ | $5.909\times 10^{-2}$
$\kappa_{Y}=0.3$, $m_{Y}=2000\text{ GeV}$ | $Ybj$ | $2b+3j$ | $1.734\times 10^{-2}$
$\kappa_{Y}=0.15$, $m_{Y}=500\text{ GeV}$ | $Ybj$ | $2b+3j$ | $3.531\times 10^{-1}$
$\kappa_{Y}=0.15$, $m_{Y}=1000\text{ GeV}$ | $Ybj$ | $2b+3j$ | $5.931\times 10^{-2}$
$\kappa_{Y}=0.15$, $m_{Y}=1500\text{ GeV}$ | $Ybj$ | $2b+3j$ | $1.513\times 10^{-2}$
$\kappa_{Y}=0.15$, $m_{Y}=2000\text{ GeV}$ | $Ybj$ | $2b+3j$ | $4.583\times 10^{-3}$
$t\bar{t}$ | $t\bar{t}$ | $2b+4j$ | $2.174\times 10^{2}$
$SingleTop$ | $tj$ | $b+3j$ | $1.346\times 10^{2}$
$tb$ | $2b+2j$ | $4.436\times 10^{0}$
$tW$ | $b+4j$ | $9.267\times 10^{-2}$
$tZj$ | $tZj$ | $b+5j$ | $1.040\times 10^{-2}$
$Dibosons$ | $WW$ | $4j$ | $2.775\times 10^{1}$
$WZ$ | $4j$ | $8.341\times 10^{0}$
$ZZ$ | $4j$ | $2.529\times 10^{0}$
$W+jets$ | $W+j$ | $3j$ | $2.540\times 10^{4}$
$W+2j$ | $4j$ | $1.062\times 10^{4}$
$W+3j$ | $5j$ | $4.329\times 10^{3}$
$W+bj$ | $b+3j$ | $4.564\times 10^{0}$
$W+bjj$ | $b+4j$ | $1.881\times 10^{0}$
$Z+jets$ | $Z+j$ | $3j$ | $6.476\times 10^{3}$
$Z+2j$ | $4j$ | $2.668\times 10^{3}$
$Z+3j$ | $5j$ | $1.086\times 10^{3}$
$Z+bj$ | $b+3j$ | $9.010\times 10^{-3}$
$Z+bjj$ | $b+4j$ | $9.828\times 10^{-3}$

The signal process $pp\rightarrow Ybj$ (which includes both $Y$ and $\bar{Y}$
VLQs, as well as $b$ and $\bar{b}$ quarks) process suffers from various
backgrounds and one of them is clearly $t\bar{t}$ production process. When
considering production and all hadronic decays, the $t\bar{t}$ process mimics
the most characteristic feature of the signal process which leads the
existence of two b-jets $(t\bar{t}\rightarrow W^{+}bW^{-}\bar{b}\to 4j+2b)$.

$t\bar{t}$ pair production is an important background in most BSM searches.
Top quark pairs generally decay into two $W$ bosons and a pair of $b$ quarks.
Then, the final state contains at least six hadronic jets. Moreover due to the
high cross section of the $t\bar{t}$ background, it shows presence at every
relevant region and resembles in many aspects of signal process, hence it must
be trimmed by applying proper event selection and analysis cuts. $t\bar{t}$
background and signal processes differ in some respects such as the existence
of more energetic jets and higher hadronic transverse energy for the signal.
However, presence of two top quarks with high cross section is obviously
problematic. As a direct consequence of this result, we expect an affection at
top mass reconstruction. Hence our cuts are decided to reduce its dominance.

Single top production also matters and needed to be handled well using proper
cuts. Since it has large cross section and some modes of that background are
very close to signal process besides due to high mass of top quark, mass
reconstruction for these processes resemble in some aspects with the signal.
Single top quark production $tq$, $tb$ and $tW$ samples, in which top quark
decays to a $W$ boson (decaying hadronically) and a $b$ quark ($t\rightarrow
Wb\rightarrow q\bar{q^{\prime}}b$), have been produced. The $tZj$ sample has
been produced with the top quark decaying hadronically ($t\rightarrow
Wb\rightarrow q\bar{q^{\prime}}b$) and $Z$ boson decaying to two quarks
($Z$$\rightarrow q\bar{q}$($q\neq b$)).

We produce $W$+$jets$ and $Z$+$jets$ samples which include only hadronic $W$
or $Z$ boson decays, and a hadronic transverse energy cut $H_{T}>400$ GeV
(where $H_{T}$ is the scalar sum of transverse momentum of selected jets) is
applied in the analysis to reduce these background.

The dibosons $WW$, $ZZ$ and $WZ$ samples have been generated with hadronic
decays leading at least four hadronic jets, which can be suppressed according
to the selections of two jets invariant mass interval.

At largely, the backgrounds mentioned after single top are less problematic.
Although some of them have large cross sections, with the invariant mass
reconstruction, cut down on the signal dominant region. Futhermore, transverse
momentum of the leading jet have significantly different behaviour for the
signal. On the whole, background can be kept under control especially at high
mass region.

The detector response to the signal and background events is simulated using a
detailed description of the HL-LHC detector card implementation in Delphes
v3.4 (key-15, ). All events are analyzed by developing an analysis macro with
Root v6 (key-16, ) . The kinematical distributions are normalized to the
number of events, which is defined to be the cross-section times luminosity
for each process under consideration.

##  VI Analysis And Results

The analysis targets $Wb$ events with $W$ boson decaying hadronically in order
to reconstruct vectorlike $Y$ quarks (full hadronic mode: $Y\to bW\to bjj$).
In this mode event selection requires at least five jets and at least one of
them being $b$-tagged, where hadronically decaying $W$ boson identified from
dijet invariant mass. Our all hadronic analysis focuses on final states at
least five (small radius) jets, one of them can be identified as forward jet
and one of the $b$-tagged jet can be considered collinear jet, other two jets
can be considered as central jets reconstructing the $W$ boson mass. These two
jets together with one leading ($b$-tagged) jet reconstruct vectorlike $Y$
quark invariant mass. The final discriminant variables are chosen as leading
jet $p_{T}$, hadronic transverse energy $H_{T}$ (We define $H_{T}$ as the
hadronic activity typically associated with single vectorlike $Y$ quark
production, where this variable is defined as the scalar sum of the transverse
momentum of selected jets in the events.), and the reconstructed mass $m_{Y}$
of the vector like $Y$quark candidate.

The normalized distributions of the number of jets in signal events (with mass
$m_{Y}=1500$ GeV and coupling $\kappa_{Y}=0.5$), and also in the relevant
backgrounds are given in Fig. 4. The signal has mostly four and five jets in
the final state, while the background $tZj$ and $t\bar{t}$ dominates in this
region, event selection with the requirement $n_{jet}\geqslant 5$ reduce other
backgrounds.

Figure 4: The normalized distribution of the number of jets (jet size) for
signal (with $m_{Y}=1500$ GeV and $\kappa_{Y}=0.5$), and also for the relevant
backgrounds in all hadronic channels.

One of the important background for our analysis is the top pair production
$t\bar{t}$ due to its high cross section. Here, we give the $p_{T}$ and $\eta$
distributions of the jets ($p_{T}$ ordered) for this process without any cuts
to reflect it’s main characteristics. Hence, we have the chance to exploit its
main differences. As we can see from a quick comparison between the signal
distributions, discriminating differences occur at the distributions of the
high $p_{T}$ jets.

The transverse momentum ($p_{T}$) and pseudo-rapidity ($\eta$) distributions
of five jets for signal, for $m_{Y}=1500$ GeV and $m_{Y}=2000$ GeV
($\kappa_{Y}=0.5$ for both) and also for top pair, associated $tZj$
production, single top, $W$+$jets$, $Z$+$jets$ and dibosons background are
presented in Fig. A.1 \- A.8, respectively.

The transverse momentum ($p_{T}$) and pseudo-rapidity ($\eta$) distributions
of jets for $W+jets$ (including $W+j,W+2j,W+3j,W+bj,W+bjj$) background are
presented in Fig. A.6. The generator level cuts are applied to these kinematic
distributions of the final state jets. Comparing the signal and the $W+jets$
background distributions, it can be seen that 1st and 2nd jets have different
distributions of the $p_{T}$ and $\eta$. A similar discussion is also valid
for the $Z+jets$ background when the hadronic mode is considered.

The transverse momentum and pseudo-rapidity distributions of all jets for
whole diboson background (including $WW,WZ,ZZ$) before performing the cuts are
shown in Fig. A.8. As it can be seen from Fig. A.8, the pseudo-rapidity
distribution of the jets are symmetrical with respect to the beam axis and the
tails of the distribution for the 5th jet are getting longer in both forward
and backward direction. The transverse momentum distribution falls rapidly for
high $p_{T}$ values.

In addition to the $p_{T}$ and $\eta$ cuts on the selected final state jets, a
transverse hadronic energy cut $H_{T}>400$ GeV is applied in the analysis,
where $H_{T}$ is the scalar sum of the $p_{T}$ of all selected jets. The
angular distance (or angular separation) between two jets defines how much two
jets are moving in the same direction, and it is usually denoted as $\Delta
R=\sqrt{\Delta\eta^{2}+\Delta\phi^{2}}$, where a cut $\Delta R(j,j)>0.4$ is
applied to the selected events. The preselection and cut flow for the analysis
are given in Table 2.

Cut efficiencies ($\%$) for different signal benchmarks ($\kappa_{Y}=0.5$ and
$m_{Y}=1000$, $1500$, $2000,$ $2500$ GeV) and backgrounds in Table 3 have been
obtained by applying the cuts explained in Table 2, where the leading jet
$p_{T}$ variable has been used.

The invariant mass distribution $m_{Y}$ for all signals and background are
presented in Fig. 5. The leading jet $p_{T}$ distributions for all signals and
backgrounds are shown in Fig. 6. The scalar sum $p_{T}$ for all signals and
backgrounds are shown in Fig. A.9. The ratio of the scalar sum $p_{T}$ to the
total scalar $H_{T}$, $H_{T}^{R}$, in the event is shown in Fig. 7. Fig. 8 and
Fig. 9 show the invariant mass distributions of the vectorlike $Y$quark for
$\kappa_{Y}=0.3$ and $\kappa_{Y}=0.15$, respectively. The distributions shown
in Fig. 5-9 and Fig. A.9 are normalized to the expected number of events,
defined as the cross-section of related process multiplied by the integrated
luminosity of $L_{int}=1000$ fb-1.

Figure 5: The invariant mass distribution $m_{Y}$ for all signals
($\kappa_{Y}=0.5$) and backgrounds. For $m_{Y}=1000$ and $1500$ GeV (upper),
for $m_{Y}=2000$ and $2500$ GeV (lower). Figure 6: The leading jet $p_{T}$
distributions for all signals ($\kappa_{Y}=0.5$) and backgrounds.  Figure 7:
The ratio of the scalar sum $p_{T}$ to the total scalar $H_{T}$.
$\kappa_{Y}=0.5$ for all signals.

Figure 8: The invariant mass distribution $m_{Y}$ for all signals
($\kappa_{Y}=0.3$) and backgrounds. For $m_{Y}=500$ and $1000$ GeV (upper),
for $m_{Y}=1500$ and $2000$ GeV (lower).

Figure 9: The invariant mass distribution $m_{Y}$ for all signals ($\kappa_{Y}=0.15$) and backgrounds. For $m_{Y}=500$ and $1000$ GeV (upper), for $m_{Y}=1500$ and $2000$ GeV (lower). Table 2: Preselection and set of cuts for the analysis of signal and background events.  Cuts | Definition
---|---
Cut-0 |  Preselection: number of jets $N_{j}\geqslant 5$
Cut-1 |  Leading jet $p_{T}>200$ GeV and $|\eta|<2.0$
Collinear jet $p_{T}>30$ GeV and $|\eta|<3.0$
Other jets $p_{T}>30$ GeV and $|\eta|<2.5$
Cut-2 |  Number of b tagged jet $N_{b}\geqslant 1$
Cut-3 |  Angular separation between two jets $\Delta R(j,j)>0.4$
Cut-4 |  Scalar sum of $p_{T}$ of jets $H_{T}>400$ GeV
Cut-5 |  Invariant mass interval of two jets $|m_{jj}-m_{W}|<20$ GeV
Cut-6 |  Invariant mass interval of three jets $|m_{bjj}-m_{Y}|<0.2m_{Y}$ GeV
Table 3: Cut efficiencies (%) for different signal benchmarks ( $\kappa_{Y}=0.5$ and $m_{Y}=1000,1500,2000,2500$ GeV) and backgrounds.  | Cut-1 | Cut-2 | Cut-3 | Cut-4 | Cut-5
---|---|---|---|---|---
$m_{Y}=1000$ GeV | $30.28$ | $25.20$ | $25.03$ | $25.03$ | $7.60$
$m_{Y}=1500$ GeV | $30.41$ | $24.42$ | $24.35$ | $24.35$ | $4.09$
$m_{Y}=2000$ GeV | $31.84$ | $25.31$ | $25.27$ | $25.27$ | $2.58$
$m_{Y}=2500$ GeV | $33.52$ | $25.71$ | $25.71$ | $25.71$ | $2.15$
$t\bar{t}$ | $6.58$ | $4.14$ | $4.14$ | $4.14$ | $2.09$
$tb$ | $7.63$ | $6.80$ | $6.71$ | $6.71$ | $3.22$
$tW$ | $13.78$ | $6.31$ | $6.28$ | $6.28$ | $3.02$
$tj$ | $2.88$ | $1.87$ | $1.87$ | $1.87$ | $1.01$
$tZj$ | $8.37$ | $4.25$ | $4.25$ | $4.25$ | $1.72$
$Wbj$ | $3.13$ | $2.30$ | $2.30$ | $2.30$ | $0.21$
$Wbjj$ | $4.48$ | $2.65$ | $2.65$ | $2.65$ | $0.71$
$Wj$ | $1.83$ | $0.61$ | $0.61$ | $0.61$ | $0.61$
$Wjj$ | $4.53$ | $0.40$ | $0.40$ | $0.40$ | $0.10$
$Wjjj$ | $6.65$ | $0.76$ | $0.76$ | $0.76$ | $0.27$
$Zbj$ | $6.67$ | $3.74$ | $3.74$ | $3.74$ | $1.72$
$Zbjj$ | $11.22$ | $5.83$ | $5.83$ | $5.83$ | $2.49$
$Zj$ | $3.29$ | $0.00$ | $0.00$ | $0.00$ | $0.00$
$Zjj$ | $4.51$ | $0.88$ | $0.88$ | $0.88$ | $0.10$
$Zjjj$ | $6.53$ | $0.69$ | $0.69$ | $0.69$ | $0.20$
$WW$ | $7.10$ | $0.81$ | $0.81$ | $0.81$ | $0.00$
$ZZ$ | $5.03$ | $0.41$ | $0.41$ | $0.41$ | $0.27$
$WZ$ | $5.54$ | $0.66$ | $0.66$ | $0.66$ | $0.66$

For calculating signal significance ($SS$) as incorporated in (key-17, ;
key-18, ; key-19, ), we use the invariant mass interval defined as
$|m_{bjj}-m_{Y}|<0.2m_{Y}$. The expected signal significance is given in terms
of signal ($S$) and background events ($B$).

In the special case of well known background, the statistical significance
($SS$)

$SS=\sqrt{2\left[(S+B)\ln\left(1+\frac{S}{B}\right)-S\right]}$

which would further reduces to $S/\sqrt{B}$ in the limit of large background
events. The signal significances for the couplings $\kappa_{Y}=0.5$,
$\kappa_{Y}=0.3$ and $\kappa_{Y}=0.15$ at the projected integrated
luminosities of $L_{int}=300$ fb-1, $1000$ fb-1 and $3000$ fb-1 are given in
Table 3, 4 and 5, respectively. These results are obtained for the invariant
mass $m_{Y}$ variable after applying the cuts as in Table 2 to all signals and
backgrounds event samples.

Table 4: Statistical significance ($SS$) for different signal benchmarks ( $\kappa_{Y}=0.5$ and $m_{Y}=1000,1500,2000,2500$ GeV) and different integrated luminosity projections at HL-LHC.  Signal |  $L_{int}($fb-1)
---|---
$300$ | $1000$ | $3000$
$\kappa_{Y}=0.5$, $m_{Y}=1000\text{ GeV}$ | $9.13$ | $16.67$ | $28.88$
$\kappa_{Y}=0.5$, $m_{Y}=1500\text{ GeV}$ | $3.73$ | $6.80$ | $11.78$
$\kappa_{Y}=0.5$, $m_{Y}=2000\text{ GeV}$ | $1.33$ | $2.43$ | $4.21$
$\kappa_{Y}=0.5$, $m_{Y}=2500\text{ GeV}$ | $0.50$ | $0.92$ | $1.59$
Table 5: The same as Table 4, but for $\kappa_{Y}=0.3$ and $m_{Y}=500,1000,1500,2000$ GeV. Signal |  $L_{int}($fb-1)
---|---
$300$ | $1000$ | $3000$
$\kappa_{Y}=0.3$, $m_{Y}=500\text{ GeV}$ | $13.54$ | $24.73$ | $42.83$
$\kappa_{Y}=0.3$, $m_{Y}=1000\text{ GeV}$ | $3.22$ | $5.88$ | $10.18$
$\kappa_{Y}=0.3$, $m_{Y}=1500\text{ GeV}$ | $0.72$ | $1.32$ | $2.28$
$\kappa_{Y}=0.3$, $m_{Y}=2000\text{ GeV}$ | $0.50$ | $0.92$ | $1.59$
Table 6: The same as Table 4, but for $\kappa_{Y}=0.15$ and $m_{Y}=500,1000,1500,2000$ GeV. Signal |  $L_{int}($fb-1)
---|---
$300$ | $1000$ | $3000$
$\kappa_{Y}=0.15$, $m_{Y}=500\text{ GeV}$ | $3.60$ | $6.58$ | $11.39$
$\kappa_{Y}=0.15$, $m_{Y}=1000\text{ GeV}$ | $0.77$ | $1.41$ | $2.45$
$\kappa_{Y}=0.15$, $m_{Y}=1500\text{ GeV}$ | $0.17$ | $0.31$ | $0.53$
$\kappa_{Y}=0.15$, $m_{Y}=2000\text{ GeV}$ | $0.15$ | $0.28$ | $0.48$

The statistical significance $SS$ as a function of mass of vectorlike $Y$
quark for $L_{int}=300$ fb-1, $1000$ fb-1 and $3000$ fb-1 are shown in Fig.
10, 11 and 12. In these figures, the coupling values $\kappa_{Y}=0.5$,
$\kappa_{Y}=0.3$ and $\kappa_{Y}=0.15$ are chosen as benchmarks, respectively.
For the integrated luminosity projection of $3000$ fb-1, lower limits for
$m_{Y}$ are obtained as $2350$ GeV at $2\sigma$ significance level, $2150$ GeV
at $3\sigma$ (observability) and $1900$ GeV at $5\sigma$ (discovery)
significances. For other couplings $\kappa_{Y}=0.3$ and $\kappa_{Y}=0.15$ the
lower limits for $m_{Y}$ are found as $1550$ GeV and $1075$ GeV at $2\sigma$
significance level at $L_{int}=3000$ fb-1, respectively.

Figure 10: The statistical significance $SS$ for $\kappa_{Y}=0.5$ as a
function of vectorlike $Y$ quark mass $m_{Y}$ for integrated luminosity
projections (solid lines, red to green) $L_{int}=300$ fb-1, $1000$ fb-1 and
$3000$ fb-1 at HL-LHC. The horizontal dashed lines correspond to $2\sigma$,
$3\sigma$ and $5\sigma$ significances (lower to upper), respectively.  Figure
11: The same as Fig. 10, but for $\kappa_{Y}=0.3$. Figure 12: The same as Fig.
10, but for $\kappa_{Y}=0.15$.

##  VII Conclusions

We have studied single production of vectorlike $Y$ quark and its subsequent
decay ($Y\to Wb$) with phenomenological interpretations in the context of well
motivated VLQ model framework. Signal significances have been obtained
depending on the parameter space (mass in the range $1000-2500$ GeV for
coupling $\kappa_{Y}=0.5$ and $500-2000$ GeV for couplings $\kappa_{Y}=0.3$
and $\kappa_{Y}=0.15$) using all hadronic mode (at least five jets final
state). We find a significant coverage of the signal parameter space and
distinguish the signal for a mass up to $2350$ GeV for an integrated
luminosity of $3000$ fb-1 at the HL-LHC. We have used the criteria for
projected discovery (key-20, ) sensitivities to the couplings and masses. For
signal significance $2\sigma$ level we find vectorlike $Y$ quark attainable
mass limits as $1800$ GeV, $2100$ GeV and $2350$ GeV at projected integrated
luminosities $L_{int}=300$ fb${}^{-1},$$1000$ fb-1 and $3000$ fb-1,
respectively. For other coupling values $\kappa_{Y}=0.3$ ($\kappa_{Y}=0.15$)
we find attainable mass limits as $1150$ GeV ($700$ GeV), $1350$ GeV ($900$
GeV) and $1550$ GeV ($1075$ GeV) at the same integrated luminosity
projections. The systematics are known to affect the discovery reaches.
However, given the difficulty to perform a precise estimation of systematics
for projected results, the significances have been computed without
considering the systematic uncertainties (an optimistic case) for the
estimation of the potential of HL-LHC. These results can also be used for
constraining the models (key-21, ; key-22, ) predicting new heavy quarks with
exotic charge as $-4/3$. This analysis shows that HL-LHC could discover a wide
range of parameter space of vectorlike $Y$ quark models.

###### Acknowledgements.

The numerical calculations reported in this paper were partially performed at
TUBITAK ULAKBIM, High Performance and Grid Computing Center (TRUBA resources).
The work of O.C. was supported in part by the Turkish Atomic Energy Authority
(TAEA) under grant No. 2020TAEK(CERN)A5.H1.F5-25.

## Appendix A

We present following Table A.1 and Figures A.1 \- A.9 to avoid detailed
information in the main text.

Table A.1: The leading order (LO) signal cross section values (in pb) for the process $pp\rightarrow Ybj+X$ at the HL-LHC with $\sqrt{s}=14$ TeV.  $m_{Y}(\text{GeV})$ | $\kappa_{Y}=0.5$ | $\kappa_{Y}=0.2$ | $\kappa_{Y}=0.1$ | $\kappa_{Y}=0.05$
---|---|---|---|---
$1000$ | $1.08\times 10^{0}$ | $1.72\times 10^{-1}$ | $4.30\times 10^{-2}$ | $1.08\times 10^{-2}$
$1100$ | $8.01\times 10^{-1}$ | $1.28\times 10^{-1}$ | $3.21\times 10^{-2}$ | $8.01\times 10^{-3}$
$1200$ | $6.04\times 10^{-1}$ | $9.63\times 10^{-2}$ | $2.42\times 10^{-2}$ | $6.04\times 10^{-3}$
$1300$ | $4.61\times 10^{-1}$ | $7.37\times 10^{-2}$ | $1.85\times 10^{-2}$ | $4.60\times 10^{-3}$
$1400$ | $3.55\times 10^{-1}$ | $5.68\times 10^{-2}$ | $1.41\times 10^{-2}$ | $3.55\times 10^{-3}$
$1500$ | $2.76\times 10^{-1}$ | $4.41\times 10^{-2}$ | $1.10\times 10^{-2}$ | $2.75\times 10^{-3}$
$1600$ | $2.15\times 10^{-1}$ | $3.44\times 10^{-2}$ | $8.56\times 10^{-3}$ | $2.14\times 10^{-3}$
$1700$ | $1.69\times 10^{-1}$ | $2.70\times 10^{-2}$ | $6.73\times 10^{-3}$ | $1.69\times 10^{-3}$
$1800$ | $1.33\times 10^{-1}$ | $2.13\times 10^{-2}$ | $5.30\times 10^{-3}$ | $1.33\times 10^{-3}$
$1900$ | $1.06\times 10^{-1}$ | $1.69\times 10^{-2}$ | $4.21\times 10^{-3}$ | $1.05\times 10^{-3}$
$2000$ | $8.44\times 10^{-2}$ | $1.34\times 10^{-2}$ | $3.37\times 10^{-3}$ | $8.41\times 10^{-4}$
$2100$ | $6.72\times 10^{-2}$ | $1.07\times 10^{-2}$ | $2.70\times 10^{-3}$ | $6.71\times 10^{-4}$
$2200$ | $5.39\times 10^{-2}$ | $8.62\times 10^{-3}$ | $2.16\times 10^{-3}$ | $5.40\times 10^{-4}$
$2300$ | $4.34\times 10^{-2}$ | $6.95\times 10^{-3}$ | $1.74\times 10^{-3}$ | $4.34\times 10^{-4}$
$2400$ | $3.49\times 10^{-2}$ | $5.62\times 10^{-3}$ | $1.40\times 10^{-3}$ | $3.50\times 10^{-4}$
$2500$ | $2.83\times 10^{-2}$ | $4.54\times 10^{-3}$ | $1.13\times 10^{-3}$ | $2.84\times 10^{-4}$
$2600$ | $2.29\times 10^{-2}$ | $3.67\times 10^{-3}$ | $9.11\times 10^{-4}$ | $2.29\times 10^{-4}$
$2700$ | $1.86\times 10^{-2}$ | $2.98\times 10^{-3}$ | $7.42\times 10^{-4}$ | $1.85\times 10^{-4}$
$2800$ | $1.51\times 10^{-2}$ | $2.42\times 10^{-3}$ | $6.02\times 10^{-4}$ | $1.51\times 10^{-4}$
$2900$ | $1.23\times 10^{-2}$ | $1.97\times 10^{-3}$ | $4.91\times 10^{-4}$ | $1.23\times 10^{-4}$
$3000$ | $1.01\times 10^{-2}$ | $1.60\times 10^{-3}$ | $4.00\times 10^{-4}$ | $9.98\times 10^{-5}$

Figure A.1: Transverse momentum (upper) and pseudo-rapidity (lower)
distribution of five jets for signal with $m_{Y}=1500$ GeV and
$\kappa_{Y}=0.5$.

Figure A.2: The same as Fig. A.1, but for signal with $m_{Y}=2000$ GeV and
$\kappa_{Y}=0.5$.

Figure A.3: Transverse momentum (upper) and pseudo-rapidity (lower)
distributions of five jets for top-pair background.

Figure A.4: Transverse momentum (upper) and pseudo-rapidity (lower)
distributions of five jets for associated tZj production background.

Figure A.5: Transverse momentum (upper) and pseudo-rapidity (lower)
distributions of five jets for single top background.

Figure A.6: Transverse momentum (upper) and pseudo-rapidity (lower)
distributions of five jets for $W+jets$ background.

Figure A.7: Transverse momentum (upper) and pseudo-rapidity (lower)
distributions of five jets for $Z+jets$ background.

Figure A.8: Transverse momentum (upper) and pseudo-rapidity (lower)
distributions of five jets for diboson background.

Figure A.9: The scalar sum $p_{T}$ for all signals ($\kappa_{Y}=0.5$) and
backgrounds. For $m_{Y}=1000$ and $1500$ GeV (upper), for $m_{Y}=2000$ and
$2500$ GeV (lower).

## References

  * (1) Apollinari G et al. 2017 High-Luminosity Large Hadron Collider (HL-LHC) Technical Design Report V. 0.1 (Geneva: CERN). https://doi.org/10.23731/CYRM-2017-004.
  * (2) M. Buchkremer, G. Cacciapaglia, A. Deandrea, L. Panizzi, Model Independent Framework for Searches of Top Partners, hep-ph/1305.4172.
  * (3) J. A. Aguilar-Saavedra, R. Benbrik, S. Heinemeyer and M. Perez-Victoria, Handbook of vectorlike quarks: Mixing and single production, Phys. Rev. D 88 (2013) 094010, arXiv: 1306.0572 [hep-ph].
  * (4) ATLAS Collaboration, Search for single production of vectorlike quarks decaying into Wb in pp collisions at $\sqrt{s}=13$ TeV with the ATLAS detector, JHEP05, 164 (2019).
  * (5) CMS Collaboration, Search for single production of vectorlike quarks decaying into a b quark and a W boson in proton–proton collisions at $\sqrt{s}=13$ TeV, Physics Letters B 772, 634-656 (2017).
  * (6) Alexandra Carvalho, Stefano Moretti, Dermot O’Brien, Luca Panizzi, Hugo Prager, Single production of vector like quarks with large width at the Large Hadron Collider, Phys.Rev. D 98, 015029 (2018).
  * (7) Aldo Deandrea, Thomas Flacke, Benjamin Fuks, Luca Panizzi and Hua-Sheng Shaoc, Single production of vector-like quarks: the effects of large width, interference and NLO corrections, JHEP08(2021)107.
  * (8) J. A. Aguilar-Saavedra, Mixing with vectorlike quarks: constraints and expectations, EPJ Web Conf. 60 (2013) 16012, arXiv: 1306.4432 [hep-ph].
  * (9) Benjamin Fuks, Hua-Sheng Shao, QCD next-to-leading-order predictions matched to parton showers for vectorlike quark models, Eur. Phys. J. C (2017) 77:135. arXiv: 1610.04622 [hep-ph]
  * (10) J. Alwall, R. Frederix, S. Frixione, V. Hirschi, F. Maltoni, O. Mattelaer, H.-S. Shao, T. Stelzer, P. Torrielli, M. Zaro, The automated computation of tree-level and next-to-leading order differential cross sections, and their matching to parton shower simulations, arXiv:1405.0301 [hep-ph].
  * (11) A. Alloul, N.D. Christensen, C. Degrande, C. Duhr, B. Fuks, FeynRules 2.0 — A complete toolbox for tree-level phenomenology, Computer Physics Communications Volume 185, Issue 8, Pages: 2250-2300, (2014).
  * (12) C. Degrande, C. Duhr, B. Fuks, D. Grellscheid, O. Mattelaer, T. Reiter, UFO – The Universal FeynRules Output, Computer Physics Communications, Volume 183, Issue 6, Pages: 1201-1214, (2012).
  * (13) Torbjörn Sjöstrand, Stefan Ask, Jesper R. Christiansen, Richard Corke, Nishita Desai, Philip Ilten, Stephen Mrenna, Stefan Prestel, Christine O. Rasmussen, Peter Z. Skand, An introduction to PYTHIA 8.2, Computer Physics Communications, Volume 191, Pages 159-177 (2015).
  * (14) R.D. Ball et al., Parton distributions with LHC data, Nucl. Phys. B 867 (2013) 244 [arXiv:1207.1303] [Inspire].
  * (15) J. de Favereau, C. Delaere, P. Demin, A. Giammanco, V. Lematre, A. Mertens and M. Selvaggi, “Delphes 3, A modular framework for fast simulation of a generic collider experiment”, arXiv:1307.6346 [hep-ex].
  * (16) R. Brun, F. Rademakers, Nuclear instruments and methods in physics research section A: accelerators, spectrometers, detectors and associated equipment. New Comput. Techn. Phys. Res. V 389, 81 (1997). https://doi.org/10.1016/S0168-9002(97)00048-X.
  * (17) G. Cowan, K. Cranmer, E. Gross and O. Vitells, Asymptotic formulae for likelihood-based tests of new physics”, Eur. Phys. J. C 71, 1554 (2011) [Eur. Phys. J. C 73, 2501 (2013)] [arXiv:1007.1727 [physics.data-an]].
  * (18) Robert D. Cousins, James T. Linnemann, Jordan Tucker, Evaluation of three methods for calculating statistical significance when incorporating a systematic uncertainty into a test of the background-only hypothesis for a Poisson process, Nuclear Instruments and Methods in Physics Research A 595 (2008) 480–501.
  * (19) G. Cowan, “Two developments in tests for discovery: use of weighted Monte Carlo events and an improved measure”, Progress on Statistical Issues in Searches”, SLAC”, June 4 - 6, 2012.
  * (20) P. N. Bhattiprolu, S.P. Martin, James D. Wells, Criteria for projected discovery and exclusion sensitivities of counting experiments, arXiv:2009.07249 [physics.data-an].
  * (21) D. Chang, W. Chang, and E. Ma, Alternative interpretation of the Fermilab Tevatron top events, Phys. Rev. D 59, 091503 (1999.
  * (22) D. Chang, W. Chang, and E. Ma, Fitting precision electroweak data with exotic heavy quarks, Phys. Rev. D 61, 037301 (2000).


Broadband THz wave generation and detection
in organic crystal PNPA at MHz repetition rates

Lukasz A. Sterczewski1,∗, Jakub Mnich1, and Jaroslaw Sotor1

1Laser & Fiber Electronics Group, Faculty of Electronics, Photonics and
Microsystems, Wroclaw University
of Science and Technology, Wyb. Wyspianskiego 27, 50-370 Wroclaw, Poland

Organic nonlinear optical (NLO) crystals have emerged as efficient room-
temperature emitters of broadband THz radiation with high electric field
strengths. Although initially confined to high-pulse-energy excitation in the
millijoule range with kHz rates, recent efforts have focused on tailoring the
physical properties of organic NLO materials for compatibility with popular
MHz-rate telecommunication wavelength lasers emitting nanojoule energy pulses.
This is motivated by the large potential of such crystals for more portable
and user-safe spectroscopic systems, additionally complemented by their room-
temperature field-sensitive THz detection capabilities. In this work, we
demonstrate the MHz-rate operation of the organic crystal PNPA
((E)-4-((4-nitrobenzylidene)amino)-$N$-phenylaniline), which was recently
discovered through data mining algorithms and reported to surpass the
conversion efficiency of rivaling NLO crystals at 1 kHz repetition rate. Using
a compact 200 mW Er:fiber laser producing 18 fs pulses at 50 MHz repetition
rate, we demonstrate the THz field generation and detection capabilities of
PNPA via performing gas phase spectroscopy in the 1.3–8.5 THz range. We obtain
a dynamic range of 40 dB over 3.6 s and 70 dB over 2 hours. Our work extends
the family of organic crystals compatible with telecommunication-wavelength
excitation using nJ pulses for spectroscopy beyond 5 THz.

00footnotetext<EMAIL_ADDRESS>
##  Introduction

Organic NLO crystals have received considerable attention in the recent decade
due to their broadband THz wave generation capabilities 1. They fill an
important niche in the source spectral coverage by complementing the region
above 5 THz, which is difficult to access using typical photoconductive
switches owing to parasitic electrical parameters or insufficiently fast
carrier dynamics 2, 3, or using quantum cascade lasers (QCLs) due to photon-
phonon coupling – the _Restrahlen_ band effect 4, 5, 6. Although a plethora of
different THz generation schemes exists, organic NLO crystals feature many
practical advantages like compatibility with popular optical excitation
wavelengths at moderate pump powers, ease of fabrication (growth), and
operation at room temperature.

Similar to inorganic crystals like gallium phosphide (GaP) or zinc telluride
(ZnTe), which have traditionally been used for THz spectroscopy 7, organic NLO
crystals offer broadband phase matching 1 reaching more than 10 THz (up to 30
THz with some dips) for crystals like DAST 8, DSTMS 9, 10, OH-1 11 or BNA 12.
The main difference, however, lies in their much higher NLO coefficients (up
to 10$\times$) and larger pump wavelength flexibility. THz generation relies
in these sources on optical rectification –- a second order nonlinear process,
which requires one to match the group index at the pump wavelength (near-
infrared) with the phase (refractive) index at THz frequencies (far-infrared).
In most cases, sub-millimeter thicknesses are preferred in organic NLO
crystals due to the limited coherence length, above which the nonlinear
conversion efficiency degrades.

Figure 1: Experimental setup and principle of THz-induced lensing detection
mechanism in PNPA a) Experimental setup constituting a THz time-domain
spectrometer. (b) THz-induced lensing when no external THz field is present,
when the field is positive (inducing a diverging lens), and when the field is
negative (inducing a converging lens). In the formula, $r_{\mathrm{EO}}$ is
the electro-optic coefficient (in pm/V), and $n$ is the refractive index at
the probe wavelength (near-infrared).

What makes organic NLO crystals of even larger relevance is their bi-
functionality – they can simultaneously serve as efficient field-sensitive
detectors and emitters. Unlike in conventional electro-optic sampling systems,
the detection process does not employ a field-induced birefringence. Instead,
the Pockels effect is used, which is a linear change of the refractive index
in the presence of an external THz $E$-field. It is responsible for the
spatial modulation of the near-infrared probe beam incident on the crystal
with a $\chi^{(2)}$ nonlinearity. This effect, known as THz-induced lensing
(TIL), can be used for THz detection as proposed by Schneider et al. 13. At
high peak fields, the dynamically induced lens focuses or de-focuses the beam
so that it can be detected by a small-aperture photodetector. Note that the
Pockels effect considered here differs significantly from the Kerr effect,
which is a third-order nonlinear process ($\chi^{(3)}$) proportional to field
intensity ($|E(t)|^{2}$) that disregards the optical phase.

From a practical perspective, a system using organic NLO crystals represents
an elegant solution to broadband THz spectroscopy in amplitude and phase. One
can use a lock-in detection scheme in tandem with sensitive near-infrared
semiconductor photodetectors like InGaAs to probe very long wavelengths with
high signal-to-noise ratios (SNR). Unfortunately, the dynamic range obtainable
with organic NLO crystals is still far away from that of photoconductive
antennas 14, 10, although for the latter obtainable at much lower frequencies
(1 THz).

To circumvent the limited conversion efficiency, phase matchability, and
eliminate parasitic absorption dips in the emission spectrum, researchers have
recently employed data mining algorithms to discover and develop new nonlinear
materials for the THz range 15. A successful result of this search is the new
organic nonlinear crystal PNPA, which offers the highest intensity THz
generation surpassing that of DAST or OH-1 16 or efficient second harmonic
generation 17 and is commercially available. Its significant advantage is the
native compatibility with telecommunication wavelengths (1550 nm), where
convenient Er:fiber lasers can be used for optical excitation.

Although PNPA holds a significant potential for broadband high-resolution THz
spectroscopy, prior demonstrations have focused only the high-field THz
generation properties with spectral coverage limited to 5 THz 18. Using mJ-
energy 100 fs laser pulses with kHz repetition rates from an optical
parametric amplifier (OPA), a nonlinear conversion efficiency exceeding 4% has
been demonstrated. Since OPA-based sources have rather large footprints and
require many safety precautions, for widespread adoption of PNPA greatly
desired is compatibility with compact femtosecond fiber lasers operating at
telecommunication wavelengths with MHz rates.

In this work, we demonstrate not only the THz generation capabilities of PNPA
in the low-power (200 mW average) 1550 nm excitation regime, but also uncover
its previously unexplored field-sensitive detection capabilities. Despite the
lowered conversion efficiency, pumping with nJ-level energy pulses with 50 MHz
repetition rate ($\sim$5$\times$104 higher than reported to date) and 18 fs
duration obtainable directly from the Er:fiber amplifier has enabled us to not
only to double the coverage to above 10 THz, but also uncover quasi-discrete
emission peaks at 18 THz and 23-24 THz. The spectral dynamic range obtainable
over a second timescale reaches 40 dB, which improves to 70 dB upon coherent
averaging over 2 hours. Our work extends the family of bi-functional organic
NLO crystals like DAST or DSTMS 10 suitable for broadband THz time-domain
spectroscopy using compact telecommunication wavelength pump sources.

##  Results and Discussion

### 2.1. Broadband THz generation

To study the emission and previously unknown detection capabilities of the
PNPA crystal using nJ energy pulses, we inserted two pieces of crystal (from
Terahertz Innovations, USA, 0.8 mm thickness) into a home-made THz time-domain
spectrometer (THz-TDS) to serve as a THz emitter and detector. For comparative
purposes, we also used a DSTMS crystal
((4-$N$,$N$-dimethylamino-4’-$N$’-methyl-stilbazolium
2,4,6-trimethylbenzenesulfonate) to detect THz waveforms produced by PNPA and
vice versa. The experimental setup is shown in Fig. 1a, where a femtosecond
laser system with a 18 fs pulse duration, 200 mW of average power, and 50 MHz
repetition rate is collimated to a beam waist of 2 mm (1/e2) and next split by
a thin 90/10 beam splitter into pump and probe beams. Figure 1b plots the
concept of TIL detection for convenience.

Figure 2: Er:fiber laser system used as a near-infrared pump for the
generation and detection of THz radiation using PNPA. (a) Optical spectrum
after the Er amplifier with characteristic features of self-phase modulation
responsible for spectral broadening. (b) Intensity autocorrelation along with
a sech2 fit indicating a pulse duration of 18 fs.

To ensure the broadest spectral coverage in the THz, it is essential to keep
the near-infrared excitation pulse as short as possible. Fig. 2a shows the
measured broadband emission spectrum owing its shape to self-phase modulation
(SPM) occurring in the amplifier. A result of in-fiber post-compression
following the SPM process is the autocorrelation shown in Fig. 2, which is
clean and almost free of side lobes signifying most of the power concentrated
in the main pulse. To excite the PNPA emitter, the pump arm is mechanically
modulated by an optical chopper at 5 kHz and delayed by a precise mechanical
stage moving a hollow roof prism mirror relative to the probe arm. To focus
the beam onto the crystal, we use an $f$=75 mm achromatic doublet lens (glass
materials with opposite-sign group velocity dispersions), which produces a
2$w_{\mathrm{opt}}$=75 $\upmu$m wide focal spot with a 5.6 mm depth of field.
The latter allows us to move the crystal relative to the first off-axis
parabolic mirror without significantly changing the pump fluence. Upon
illuminating the PNPA crystal, the pump beam is blocked by a 5 mm thick
germanium window serving as a long-pass filter (LPF). The THz pulse passes
through the filter to be next imaged onto the detector crystal using 4 gold-
coated off-axis parabolic mirrors ($f$=50.8 mm). The last mirror preceding the
detection stage has a hole allowing to spatially overlap the THz pulse with
the probe.

### 2.2. Broadband THz detection

For the detection, we weakly focus the probe beam to $\sim$200 $\upmu$m by a
thin, $f$=200 mm plano-convex lens, which provides a quasi-collimated beam
with a $\sim$40 mm depth of field. Empirically, 50 mm away from the crystal we
found the beam diameter to increase its diameter from the focus by $\sqrt{2}$
to a beam waist of 300 $\upmu$m, which defines our empirical Raleigh range. We
placed there an unbiased InGaAs detector with a 100 $\upmu$m pinhole to
measure only the central part of the beam intensity, as emphasized in the
seminal work in Ref. 13. Such a simple detection geometry avoids the need for
sophisticated beam telescopes to lower the beam diameter and has proven to
provide quality THz spectroscopy data in agreement the HITRAN 19 spectroscopic
database (discussed later in the text). For convenience, we provide a
conceptual drawing of TIL-based detection scheme in Fig. 1b. Since the AC
component of the detected photocurrent lies in the nA range, we used a low-
noise transimpedance amplifier (DHPCA-100) to convert the detected current
into voltage detected in a homodyne configuration using a lock-in amplifier.

Figure 3: THz spectroscopy experiment wherein the measurement chamber was
purged with argon. (a) Time-domain waveform normalized to the water-vapor free
THz signal. (b) Corresponding THz power spectra offset by 20 dB for clarity.
The propagation environment was technical argon and laboratory air (15%
relative humidity), 0.43% H2O by volume.

### 2.3. THz spectroscopy of humid air

To test the spectroscopic capabilities of PNPA, we enclosed the time-domain
spectrometer in a hermetic chamber equipped with a humidity sensor. First, we
recorded the THz spectrum without purging the propagation environment followed
by its intensive purging with argon. When water vapor lines almost completely
disappeared from the THz spectrum, suggesting a relative humidity lower than
2%, we acquired a second THz spectrum. Results of the experiment are plotted
in Fig. 3. The acquisition delay range was 30 ps, which corresponds to a
theoretical optical resolution of 33.3 GHz ($\sim$1.11 cm-1). To obtain the
spectrum, 10 step scans with an integration time of 30 ms per step were
averaged.

It is evident from Fig. 3a that removing water vapor from the measurement
chamber through purging greatly suppresses THz field oscillations after the
main pulse. It is expected since sharp absorption dips in the frequency
spectrum must correspond to long-lasting sine-wave-like oscillations in the
time domain. It should be noted that it is only the spectral region up to 10
THz (30 $\upmu$m of wavelength) that is attenuated by H20\. Weak emission
peaks around 18 THz and higher frequencies do not change in the presence of
argon instead of laboratory air.

Figure 4: Measured absorbance spectrum of humid air (15% relative humidity) in
the measurement chamber using THz-TDS plotted with a HITRAN2020 model 19.
Figure 5: Spectrogram of the PNPA emission waveform after propagation in
ambient air. (a) Sampled $E$-field. (b) Corresponding spectrogram displaying a
chirped nature of the signal – lower frequencies appear first with a gradual
broadening towards 20–30 THz (15–10 $\upmu$m of wavelength). After 1 ps, only
oscillations below 10 THz persist, which arise due to the rotational motion of
water vapor molecules along the propagation path of the THz pulse.

By subtracting base-10 logarithmic spectra from Fig. 3b, we generated the
absorbance spectrum plotted in Fig. 4. Good agreement with the HITRAN database
is obtained (USA model, mean latitude, winter. 0.3 cm-1 apodized resolution,
Voigt profile with a triangular instrumental function). This unequivocally
proves the PNPA’s ability to serve as a material for broadband THz
spectroscopy with moderate excitation powers obtainable from Er:fiber laser
systems. Minor discrepancies between the measured spectrum and the model
relate to the finite dynamic range of the instrument varying across the
spectrum (the highest at $\sim$5 THz), and insufficiently precise sample
triggering mechanism. Here we relied on the built-in encoder of the delay
line, but further extension of the spectrum to the mid-infrared range may
require interferometric, laser-based sample triggering.

### 2.4. Spectro-temporal characteristics

Access to the THz $E$-field (its scaled representation convolved with the
detector- and instrumental response function) allows us to study the spectro-
temporal characteristics of the generated signal. An analysis performed using
a short-time Fourier transform is shown in Fig. 5, where the top panel plots
the THz waveform (Fig. 5a), while the bottom includes is its spectrogram
representation Fig. 5b). As evident from the figure, the waveform is chirped
with lower frequencies reaching the detector first, followed by a gradually
broadened spectral response with frequency components above 20 THz. After
$\sim$1 ps they disappear from the spectrum and only oscillations from water
vapor molecules persist. The intensive part of the pulse is approximately
three optical cycles long and lasts less than 1 ps. The origin of the waveform
chirp can be attributed to dispersive properties of the crystal (also for the
probe) and of the Ge filter with the former being more likely to dominate. The
waveform’s temporal structure may also originate from the molecular dynamics
of the crystal, where the oscillatory part following the main peak stems from
collective vibrations of the molecular system 9. They spectrally interfere
with the broadband single-cycle THz pulse resulting directly from pump
excitation and give rise to a modulated spectrum with prominent peaks and
regions with a lack of THz power due to destructive interference 20.

It is important to note that the peak emission frequency is also dependent on
the beam waist ($w_{\mathrm{opt}}$) stimulating the emission crystal 1. For
tightly-focused beams, lower generation frequencies suffer more from
diffraction losses. For instance, for a focal spot of 75 $\upmu$m, at 2 THz we
anticipate 20 dB of THz power losses, which drops to 15 dB at 5 THz and 11 dB
at 10 THz. Doubling the beam diameter to 150 $\upmu$m lowers these values to
12 dB, 7 dB, and 4 dB, respectively. On the other hand, the peak intensity of
the pump drops quadratically with the beam diameter $w_{\mathrm{opt}}$. Since
the THz field is proportional to
$E_{\mathrm{p}}/(\frac{1}{2}\mathrm{\pi}w_{\mathrm{opt}}^{2})$, where
$E_{\mathrm{p}}$ is the pulse energy, increasing the beam diameter $k$ times
lowers the obtainable THz power by a factor or $k^{4}$. By simply doubling $k$
while keeping other pump parameters unchanged, the THz power lowers by 12 dB.
This explains why in the mJ pulse energy regime, THz emission spectra show
generally more low-frequency content 16, 18 than with nJ pulse excitation
reported here and in Ref. 10, 21. It relates to a trade-off between the beam
size limited by the damage threshold (which should be as small as possible for
high conversion efficiency) and diffraction losses affecting predominantly
lower frequencies for tightly focused beams.

Figure 6: THz spectroscopy using different generation (Tx) and detection
crystal (Rx) pairs. (a) Time-domain THz signal. (b) Power spectra plotted with
a 20 dB offset. Although DSTMS offers the broadest detection spectrum, in the
relevant spectral region (2–8 THz), PNPA can serve equally well for this
purpose.

### 2.5. Detection capabilities

Until now, we used a well-characterized, broadband detection crystal – DSTMS
to study the broadband emission properties of PNPA with sub-20-fs pulses. This
is because unlike for electrooptic sampling, retrieving the $E$-field strength
using TIL is cumbersome as the measured quantity is only the relative change
of intensity ($\Delta I/I$). Therefore, a reasonable approach to assess the
performance of organic NLO detectors is comparing new materials with well-
documented molecular standards or other inorganic crystals 13.

Following this reasoning, first we performed a series of measurements with
PNPA and DSTMS interchangeably serving as a source and detector to finally
compare PNPA with PNPA. Results of time-domain measurements are plotted in
Fig. 6a along with corresponding power spectra in Fig. 6b.

Since the time-domain waveforms are acquired over a relatively long delay
range of 20 ps, water vapor absorption features are prominent. Even though the
spectra are plotted with a 20 dB offset, emission from the PNPA crystal
detected by DSTMS offers the broadest spectral response with frequency
components filling the entire spectral region of interest (20 THz). On the
other hand, within the region of the highest spectral power (below 8 THz),
minimal differences between the spectra exists. This proves the PNPA’s
suitability for broadband THz spectroscopy.

Figure 7: Sensitivity of the detection system to axial shift when the PNPA
detection crystal was moved in 0.125 mm steps from 0 to 2.75 mm from the focus
of the off-axis parabolic mirror. (a) Stack of time-domain traces. (b) RMS
field strength calculated from (a). (c) 3 traces at arbitrary positions: 0 mm
(in focus), offset by 1 mm, and 2 mm. (d) Corresponding power spectra.

### 2.6. Conversion efficiency

To assess the performance more quantitatively, we measured the average THz
power using a broadband LiTaO3 pyroelectric detector with a HDPE window (Laser
Components, Germany). From DSTMS we measure 430 nW (after the Ge window),
which corresponds to a power conversion efficiency of
$\eta_{\mathrm{D}}$=2.4$\cdot$10-6. For PNPA we measure an average power of
170 nW, which yields $\eta_{\mathrm{P}}$=9.4$\cdot$10-7. This discrepancy, or
at least the superiority of DSTMS over PNPA in terms of THz power conversion
efficiency, can potentially be explained by the lack of broadband phase
matching in the case of PNPA with frequency components extending over 30 THz
rather than 10 THz 1. Also, prior data on PNPA as an emitter focused on the
lower-frequency THz region limited to 5 THz with most of the power
concentrated at 1.5–2.5 THz.

Pumping with mJ-energy pulses often leads to percent-level nonlinear
conversion efficiency, while for telecom-wavelength Er:fiber systems $\eta$ on
the order of 10-6 is naturally expected. We can explain it using the linear
dependence of the THz field with respect to pulse energy, which translates
into a quadratic dependence for the THz power. Accounting for the large-area
illumination ($w_{\mathrm{opt}}$=3.8 mm, $k$=50 larger than in our setup), and
0.6 mJ pulse energy (vs 3.6 nJ here, 1.7$\cdot$105 higher) we anticipate
$\sim$105 THz more power than in our setup, which is close to the mW average
THz power range, when percentage level conversion efficiency was obtained. We
want to emphasize that it is only a coarse approximation ignoring saturation
or thermal effects, yet it explains the obtained performance metrics.

### 2.7. Detection sensitivity

To characterize the detection sensitivity we measure the relative intensity
(photocurrent) change induced by the THz field produced by DSTMS and PNPA.
With average measured photocurrents (DC) lying in the $\upmu$A range, we
measure peak THz-induced photocurrent changes in the nA range. Specifically,
$\Delta I/I$=6$\cdot$10-4 for both DSTMS and PNPA serving as detection
crystals. When PNPA crystal was used for the generation and detection (with
2.5$\times$ less average power), $\Delta I/I$=3.4$\cdot$10-4. Considering the
experimental uncertainty due to alignment and nearly identical refractive
indices of DSTMS and PNPA 17, we conclude similar detection performance of
PNPA to that of DSTMS. A more accurate assessment of the detection
capabilities (including the electro-optic coefficient $r_{\mathrm{EO}}$) will
require refractometry of PNPA in tandem with calibration of the $E$-field
strength, as required by the Pockels effect formula from Fig. 1b. From the
current measurements, the most obvious observation is the lack of ultra-
broadband detection ($>$10 THz) capabilities of PNPA at telecommuncation
wavelengths. However, phase matching may be present closer to the mid-infrared
pump region, where very few organic NLO materials have been proven to operate
22.

Another practical aspect of the TIL scheme is the sensitivity to axial
position of the detection crystal with respect to the THz-focusing off-axis
parabolic mirror. In steps of 0.125 $\upmu$m, we moved the position of PNPA
relative to the point of maximum detected THz amplitude over 2.75 mm. Recorded
THz signals (Fig. 7a,c) along with the normalized RMS field strength (Fig. 7b)
and power spectra (Fig. 7d) indicate a rapid drop of the detected field
strength when the crystal was moved by more than 0.5 mm away from the focal
spot. On the other hand, even at large distances (2 mm), the THz waveform was
still detectable, albeit with $\sim$14 dB loss (0.2 of the maximum field
strength). It is important to note that almost all THz frequencies attenuate
by the same amount with little-to-no visual difference between them, which is
not expected for broadband polychromatic light detected using the TIL
technique.

The behavior observed in Fig. 7 can be explained considering the spot size of
the THz beam $w_{0}=w_{\mathrm{opt}}$/$\sqrt{2}$ (the $\sqrt{2}$ factor is
because the THz field is directly proportional to the pump intensity $I$,
hence the THz power is proportional to $I^{2}$ 10). Empirically, we measured
the THz spot size using a pyroelectric detector with a calibrated pinhole,
which in our case amounted to $\sim$150 $\upmu$m ($\pm$50 $\upmu$m). Such a
beam has a diameter for the $E$-field of $\sim$200 $\upmu$m, which doubles at
0.5 mm away from the focus, and increases its size 7 times 2 mm away
considering Gaussian beam propagation. This is in good agreement with the data
shown in Fig. 7b, but also proves that our THz imaging system (using 4 off-
axis parabolic mirrors) is not far from diffraction-limited performance if one
assumes a peak emission wavelength of 70 $\upmu$m and a collimated THz beam
diameter of $\sim$35 mm with a theoretical spot size of 130 $\upmu$m.

### 2.8. Dynamic range capabilities

Finally, we characterized the dynamic range capabilities of the THz-TDS system
employing PNPA crystals. 20 ps long traces were acquired in rapid acquisition
mode with a lock-in integration time of 1 ms. Each trace consisted of 1600
points taken every 3.6 seconds. Via cross-correlation-based alignment of THz
traces (to correct for the random jitter of the acquisition trigger event), we
obtained 2000 signals averaged in the time domain. Fig. 8 shows the dynamic
range capabilities, where a single 3.6 s-long scan provides $\sim$40 dB of
dynamic range, which improves to $\sim$70 dB upon integration. We want to
emphasize, however, that scans on a minute timescale already provide dozens of
dB of dynamic range improvement from 40 to almost 60 dB, which should already
be sufficient for many spectroscopic investigations.

Figure 8: Dynamic range capabilities of the PNPA crystal serving as a source
and detector.

##  Conclusion

In summary, we have demonstrated the nJ pulse energy compatibility of the
organic NLO crystal PNPA with an Er:fiber laser system to generate and detect
field-resolved THz waveforms. Using 18-fs long pulses with 200 mW of average
power centered at $\sim$1.56 $\upmu$m of wavelength, we have obtained 170 nW
of THz power spread over more than 20 THz with the most intensive features in
the 2–8 THz range. The emitted waveform exhibits slight chirp potentially
attributable to optical dispersion and molecular system dynamics. We have also
shown the THz-induced lensing capabilities of PNPA similar to DSTMS rendering
a peak spectral dynamic range of 40 dB over seconds and 70 dB over 2 h of
averaging. Although PNPA lacks the ultra-broadband emission characteristics of
rivaling organic NLO crystals when pumped by 1550 nm sources, it is likely to
be preferred for sub-10-THz frequency applications or to complement dips in
the emission spectra from crystals like DSTMS or DAST via THz waveform
synthesis 16.

We also want to note that ultra-broadband (70 THz coverage) Ge-based
photoconductive antennas 23 have recently been made compatible with Er:fiber
laser technology 20. As a competing platform, they show a lot of promise for
gap-less THz emission. Unfortunately, when pumped at 1550 nm (instead of the
optimal 1100 nm), they lack power between 3–10 THz. This is where PNPA can
potentially serve to fill the spectral gap.

##  Experimental Section

Calculation of Fourier spectra: Each time-domain waveform was windowed using a
Tukey window with a shape parameter of 0.1 after mean value removal. 32-times
zero padding has been used to obtain interpolated, smooth power spectra using
the fast Fourier transform (fft) algorithm.

## Acknowledgements

Acknowledgment: L. A. Sterczewski and J. Mnich acknowledge funding from the
European Union (ERC Starting Grant, TeraERC, 101117433). Views and opinions
expressed are however those of the authors only and do not necessarily reflect
those of the European Union or the European Research Council Executive Agency.
Neither the European Union nor the granting authority can be held responsible
for them. This work is supported by the use of National Laboratory for
Photonics and Quantum Technologies (NPLQT) infrastructure, which is financed
by the European Funds under the Smart Growth Operational Programme.

## Author contributions

L.A.S. conceived the idea. J.S. fabricated the femtosecond laser. L.A.S.
carried out the optical and electrical measurements with help from J.M., and
analyzed the data. L.A.S wrote the manuscript with input from J.M. and J.S..
L.A.S. acquired funding and coordinated the project.

## Conflict of interest

The authors declare no conflict of interest.

## Data availability statement

Data are available in a public, open access repository. The data that support
the findings of this study are openly available in figshare:
https://dx.doi.org/10.6084/m9.figshare.26304973.

## References

  * 1 Jazbinsek, M., Puc, U., Abina, A. & Zidansek, A.  Organic Crystals for THz Photonics.  _Applied Sciences_ 9, 882 (2019).
  * 2 Göbel, T. _et al._ Telecom technology based continuous wave terahertz photomixing system with 105 decibel signal-to-noise ratio and 35 terahertz bandwidth.  _Optics Letters_ 38, 4197 (2013).
  * 3 Turan, D. _et al._ Wavelength conversion through plasmon-coupled surface states.  _Nature Communications_ 12, 4641 (2021).
  * 4 Köhler, R. _et al._ Terahertz semiconductor-heterostructure laser.  _Nature_ 417, 156–159 (2002).
  * 5 Williams, B. S.  Terahertz quantum-cascade lasers.  _Nature Photonics_ 1, 517–525 (2007).
  * 6 Shahili, M. _et al._ Continuous-wave GaAs/AlGaAs quantum cascade laser at 5.7 THz.  _Nanophotonics_ 13, 1735–1743 (2024).
  * 7 Ferguson, B. & Zhang, X.-C.  Materials for terahertz science and technology.  _Nature Materials_ 1, 26–33 (2002).
  * 8 He, Y. _et al._ High-energy and ultra-wideband tunable terahertz source with DAST crystal via difference frequency generation.  _Applied Physics B_ 124, 16 (2018).
  * 9 Somma, C. _et al._ Ultra-broadband terahertz pulses generated in the organic crystal DSTMS.  _Optics Letters_ 40, 3404 (2015).
  * 10 Puc, U., Bach, T., Günter, P., Zgonik, M. & Jazbinsek, M.  Ultra-Broadband and High-Dynamic-Range THz Time-Domain Spectroscopy System Based on Organic Crystal Emitter and Detector in Transmission and Reflection Geometry.  _Advanced Photonics Research_ 2, 2000098 (2021).
  * 11 Vicario, C., Ruchert, C. & Hauri, C.  High field broadband THz generation in organic materials.  _Journal of Modern Optics_ 62, 1480–1485 (2015).
  * 12 Zhao, H. _et al._ Efficient broadband terahertz generation from organic crystal BNA using near infrared pump.  _Applied Physics Letters_ 5 (2019).
  * 13 Schneider, A., Biaggio, I. & Günter, P.  Terahertz-induced lensing and its use for the detection of terahertz pulses in a birefringent crystal.  _Applied Physics Letters_ 84, 2229–2231 (2004).
  * 14 Vieweg, N. _et al._ Terahertz-time domain spectrometer with 90 dB peak dynamic range.  _Journal of Infrared, Millimeter, and Terahertz Waves_ 35, 823–832 (2014).
  * 15 Valdivia-Berroeta, G. A. _et al._ Data mining for terahertz generation crystals.  _Advanced Materials_ 34, 2107900 (2022).
  * 16 Rader, C. _et al._ Custom terahertz waveforms using complementary organic nonlinear optical crystals.  _Optics Letters_ 47, 5985 (2022).
  * 17 Holland, K. M. _et al._ Characterization of organic crystals for second-harmonic generation.  _Optics Letters_ 48, 5855–5858 (2023).
  * 18 Rader, C. _et al._ A New Standard in High-Field Terahertz Generation: The Organic Nonlinear Optical Crystal PNPA.  _ACS Photonics_ 9, 3720–3726 (2022).
  * 19 Gordon, I. E. _et al._ The hitran2020 molecular spectroscopic database.  _Journal of quantitative spectroscopy and radiative transfer_ 277, 107949 (2022).
  * 20 Singh, A. _et al._ Up to 70 thz bandwidth from an implanted ge photoconductive antenna excited by a femtosecond er: fibre laser.  _Light: Science & Applications_ 9, 30 (2020).
  * 21 Puc, U., Yang, J.-A., Kim, D., Kwon, O.-P. & Jazbinsek, M.  Broadband thz wave generation in organic benzothiazolium crystals at mhz repetition rates.  _Optical Materials Express_ 13, 53–66 (2023).
  * 22 Gollner, C. _et al._ Highly efficient thz generation by optical rectification of mid-ir pulses in dast.  _APL Photonics_ 6 (2021).
  * 23 Singh, A., Pashkin, A., Winnerl, S., Helm, M. & Schneider, H.  Gapless broadband terahertz emission from a germanium photoconductive emitter.  _ACS Photonics_ 5, 2718–2723 (2018).


# LT@Helsinki at SemEval-2020 Task 12:
Multilingual or language-specific BERT?

Marc Pàmies
Emily Öhman
University of Helsinki
<EMAIL_ADDRESS>Kaisla Kajava
Jörg Tiedemann

(May 31, 2020)

###### Abstract

This paper presents the different models submitted by the LT@Helsinki team for
the SemEval 2020 Shared Task 12. Our team participated in sub-tasks A and C;
titled offensive language identification and offense target identification,
respectively. In both cases we used the so-called Bidirectional Encoder
Representation from Transformer (BERT), a model pre-trained by Google and
fine-tuned by us on the OLID and SOLID datasets. The results show that
offensive tweet classification is one of several language-based tasks where
BERT can achieve state-of-the-art results.

##  1 Introduction

00footnotetext:  This work is licensed under a Creative Commons Attribution
4.0 International License. License details:
http://creativecommons.org/licenses/by/4.0/.

The number of social media users has reached 3.5 billion, and an average of
6,000 tweets are generated every second [Mathew et al., 2019] . With such a
large volume of tweets, it seems inevitable that some use offensive language.
A study from 2014 found that 67% of social media users had been exposed to
online hate, and 21% had been the target of online hate [Oksanen et al.,
2014].

The usual approach on social media sites is to forbid hate speech in the terms
of service and to censor any inappropriate content detected by their
algorithms, but still, companies like Facebook or Twitter have been harshly
criticized for not doing enough [Del Vigna et al., 2017]. This criticism has
forced companies to try and find accurate and scalable solutions that solve
the problem of offensive language detection using automated methods. Workshops
dealing with offensive language, such as TRAC [Kumar et al., 2018b], TA-COS
[Lefever et al., 2018], or ALW1 [Waseem et al., 2017], as well as shared tasks
like GermEval [Wiegand et al., 2018], TRAC-1 [Kumar et al., 2018a], and
OffensEval 2019 [Zampieri et al., 2019b], are becoming more and more
prevalent.

The task we address here, SemEval 2020 task 12 [Zampieri et al., 2020], is
titled Multilingual Offensive Language Identification in Social Media and is
divided into the following sub-tasks:

> A.
>
>
> Offensive Language Identification in several languages (Arabic, Danish,
> English, Greek, Turkish): whether a tweet is offensive or not.
>
> B.
>
>
> Categorization of Offense Types: whether an offensive tweet is targeted or
> untargeted.
>
> C.
>
>
> Offense Target Identification: whether a targeted offensive tweet is
> directed towards an individual, a group or otherwise.

In this paper the system created by the LT@Helsinki team for sub-tasks A and C
will be described. In sub-task A we participated in all the language tracks.
For sub-task C the only language available was English. We qualified as second
in sub-task C and our submission for sub-task A ranked first for Danish,
seventh for Greek, eighteenth for Turkish, and forty-sixth for Arabic. In all
submissions we used BERT-Base models [Devlin et al., 2019] fine-tuned on each
dataset provided by the task organizers. We also experimented with random
forest with TF-IDF and other kinds of features, but the results on the
development set were not as good as with transfer learning techniques based on
pre-trained language models. We discovered that, at least for this data, the
language-specific model worked better than the multilingual.

##  2 Background

All the data used to train our classifiers was provided by the OffensEval
organizers [Zampieri et al., 2019a]. The tweets were retrieved from the
Twitter Search API and manually labeled by at least two human annotators. As a
pre-processing step, they desensitized all tweets replacing usernames and
website URLs by general tokens. The following datasets were used for the
offensive language identification problem:

Language | Training | Test | Total
---|---|---|---
Arabic | 8,000 | 2,000 | 10,000
Danish | 2,960 | 329 | 3,289
Greek | 8,743 | 1,544 | 10,287
Turkish | 31,756 | 3,528 | 35,284
Table 1: Sub-task A datasets

Some of the algorithms that can be found in the literature are random forest
[Burnap and Williams, 2015], logistic regression [Davidson et al., 2017], and
Support Vector Machine [Malmasi and Zampieri, 2018], as well as deep learning
approaches like Convolutional Neural Networks [Gambäck and Sikdar, 2017] or
Convolutional-GRU [Zhang et al., 2018]. However, in 2018 deep pre-trained
language models obtained state-of-the-art results in several NLP downstream
tasks, text classification being one of them. In particular, Google’s
Bidirectional Encoder Representations from Transformers (BERT) [Devlin et al.,
2019] stood above the rest for being deeply bidirectional and using the novel
self-attention layers from the transformer model [Vaswani et al., 2017], which
allows it to better understand a word’s context by looking at both its left
and right neighbours. BERT provides out-of-the-box pre-trained monolingual and
multilingual models that, after massive training on general corpora, can be
fine-tuned with a small amount of task-specific data and still offer excellent
performance. The results published in OffensEval’s previous edition [Zampieri
et al., 2019b] proved that BERT is well suited for the offensive language
detection task, since it was the chosen method of most of the top teams,
including the winners of sub-tasks A [Liu et al., 2019] and C [Nikolov and
Radivchev, 2019].

For the Danish dataset [Sigurbergsson and Derczynski, 2019] we used Nordic
BERT111https://github.com/botxo/nordic_bert, which is pre-trained on Danish
Wikipedia texts, Danish text from Common Crawl, Danish OpenSubtitles, and text
from popular Danish online forums. All in all, the training corpus consists of
over 90M sentences and almost 20M unique tokens. For the other languages we
used the standard BERT-Base models [Devlin et al., 2019] with no further pre-
training and only the provided datasets for each language, i.e. Arabic
[Mubarak et al., 2020], Greek [Pitenis et al., 2020], and Turkish [Çöltekin,
2020].

In sub-task C, since the language was English it was possible to use the
Offensive Language Identification Dataset (OLID), provided by the organizers
of OffensEval 2019. Despite consisting of different sub-tasks, all of them
shared the same dataset that was annotated according to a three-level
hierarchical model, so that each sub-task could use as dataset a subset of the
previous sub-task’s dataset. First, all tweets were labeled as either
offensive (OFF) or not offensive (NOT). Then, for sub-task B, all the
offensive tweets were labeled as targeted (TIN) or untargeted insults (UNT).
And finally, for the last sub-task, the third level of the hierarchy labeled
targeted insults based on who was the recipient of the offense: an individual
(IND), a group (GRP) or a different kind of entity (OTH). To illustrate this,
Figure 2 displays OLID’s label distribution.

A | B | C | Training | Test | Total
---|---|---|---|---|---
OFF | TIN | IND | 2,407 | 100 | 2,507
OFF | TIN | GRP | 1,074 | 78 | 1,152
OFF | TIN | OTH | 395 | 35 | 430
OFF | UNT | - | 524 | 27 | 551
NOT | - | - | 8,840 | 620 | 9,460
All |  |  | 13,240 | 860 | 14,100
Table 2: Distribution of label combinations in OLID.

The corpus from sub-task C in OffensEval 2019 contains 4,089 English tweets,
of which 3,876 originally belonged to the training set and the remaining 213
to the test set (Figure 2). We also had at our disposal over 9M English tweets
provided by the OffensEval 2020 organizers [Rosenthal et al., 2020]. However,
these were processed by unsupervised learning methods instead of human
annotators, which is why on this occasion each tweet was not associated with a
label but given two values: (1) the confidence that it belongs to a specific
class and (2) its standard deviation.

##  3 System Overview and Experimental Setup

In this section we discuss the system we created and setup we used for subtask
A and subtask C.

###  3.1 Sub-task A

For the binary classification problem (sub-task A) two baseline methods were
implemented and evaluated on the training Danish dataset: random forest and
BERT.

For the random forest implementation, the pre-processing steps were lower-
casing all characters, removing irrelevant punctuation marks, reducing the
length of characters that appear more than two consecutive times, and
converting hashtags to sentences by adding white spaces before every capital
letter. Tokenization was done with the TweetTokenizer tool from NLTK. The same
library was also used to perform stopword removal and word stemming. Emojis
were removed after storing their “sentiment score” as a feature employing the
emosent Python utility package [Novak et al., 2015]. We also used surface-
level features such as the number of URL tokens and @USER mentions, the total
number of characters, punctuation marks and words in each post, average word
length, percentage of capital letters, and the number of abusive terms. A
ratio of 10:1 was applied when splitting the dataset into training and
validation sets. Since the Danish dataset was relatively small, 10-fold cross-
validation was done to obtain reliable results. Otherwise the F1 scores relied
too much on which samples would fall in the validation set. The random forest
implementation from scikit-learn [Pedregosa et al., 2011] was used. The
optimal parameters were found with grid search.

System | macro-F1
---|---
All NOT | 0.465
Random Forest with TFIDF | 0.773
Multilingual BERT-Base | 0.768
Nordic BERT | 0.804
Table 3: Development results sub-task A, Danish dataset

The other approach was to apply pre-trained BERT models. After experimenting
with both the original Base version (12-layer, 768-hidden, 12-heads, 110M
parameters) and the publicly available Nordic BERT
222https://github.com/botxo/nordic_bert it was clear that the second one was
better suited for the task. This shows that further pre-training of BERT can
significantly boost performance, especially in cases like this where there is
very little data (the Danish dataset was by far the smallest of the shared
task). Our final submission was generated by the Danish version of Nordic BERT
fine-tuned on the OffensEval 2020 data, using a batch size of 32 for training
and 16 for both validation and testing. The learning rate of Adam optimizer
was set to 2e-5 and the model was trained for 4 epochs. The sequence length
was set to 128 because, even though some instances are very long (some of them
are not tweets but Reddit comments), an analysis of the length distribution
showed that only 3.4% of the training examples reached the limit after being
tokenized by BertTokenizer.

Due to time constraints we were not able to experiment in depth with all
languages available for sub-task A, which is why Danish was our main focus.
However, seeing that BERT could perform so well with almost no pre-processing,
we decided to generate results for the other three languages as well. For
Turkish and Arabic, we used the BERT-base-multilingual-cased model with a
maximum sequence length of 128, a batch size of 32, and a learning rate of
2e-5. On the other hand, for Greek we used the BERT-base-uncased model after
lowercasing and translating the entire dataset into English. In all cases we
trained the model for 4 epochs, used BertTokenizer to tokenize the tweets, and
padded and truncated the sequences to make sure each data instance had the
same length.

###  3.2 Sub-task C

The target identification problem (sub-task C, English) had an additional
level of difficulty with respect to the first task because the dataset was
highly imbalanced and composed of three classes instead of two. Thus, in this
case we focused our efforts on balancing the given dataset to prevent having a
high number of misclassifications from instances of the minority class (which
would notably affect the resulting macro-F1 score). In order to overcome the
class imbalance problem, we trained our model on all the data from 2019 and
some additional instances from the non-majority classes (GRP and OTH) of the
2020 dataset. Only the 300 OTH instances and 237 GRP instances of highest
confidence were added in order to slightly increase the balance of the
dataset. We experimented with different thresholds to select more samples but
at the end we decided to keep the value low to ensure that all tweets used for
training are tagged correctly. Finally, to overcome the class imbalance
problem, an over-sampling technique with replacing was applied. We simply
produced copies of instances from the minority classes to end up with a
totally balanced dataset of 11,628 tweets (3,876 for each class). Under-
sampling was not an option because then we would be facing a data scarcity
problem, and we experimented with ratios other than 1:1:1 obtaining promising
results but not significantly better than the proposed approach. Another
interesting approach, which was chosen by the winners of last year’s edition
[Nikolov and Radivchev, 2019], is to modify the classification thresholds
(i.e. lower the thresholds for classes OTH and GRP) to get less new examples
classified as the majority class.

Then, the now balanced dataset was used as input for a BERT-base-uncased model
with a maximum sequence length of 128 and batch sizes of 32, 16 and 8 for
training, validation and prediction respectively. The learning rate was set to
2e-5 and the training lasted 4 epochs. In this case we did not perform any
pre-processing step other than lowercasing since none of the attempts (i.e.
translating emojis to sentences, splitting hashtags into separate words etc.)
significantly boosted performance. Moreover, we believe that for this specific
task it might not be a good idea to remove @USER mentions or certain stop
words since they can carry valuable information for target identification.

###  3.3 Unfruitful Attempts at Performance Improvement

| pos | neg | ang | ant | dis | fea | joy | sad | sur | tru
---|---|---|---|---|---|---|---|---|---|---
NOT | 51.476 | 79.136 | 30.028 | 17.623 | 29.685 | 47.055 | 16.588 | 25.843 | 10.753 | 26.785
OFF | 68.640 | 110.350 | 41.481 | 23.474 | 37.543 | 58.420 | 22.374 | 46.046 | 9.190 | 26.838
Table 4: Emotion word distribution in offensive and non-offensive messages.

Offensive tweets and posts had statistically significant amounts of emotion-
laden words, including positive ones. We arrived at these results by
augmenting the data (see Ohman [Öhman et al., 2016]) using the NRC Emotion
Lexicon [Mohammad and Turney, 2013] for the languages in question and then
comparing the normalized word counts per offensive post. We did not take
context or negation into account. Unfortunately, we were unable to
meaningfully utilize this information to improve the performance of our
system. Even though offensive tweets were more likely to include emotion words
and more of them, a significant number of offensive tweets contained no
emotion words and many non-offensive tweets did contain them. Nonetheless, the
results were quite interesting so we wanted to share them here in the hopes
that they might be useful to someone else, perhaps for the OffensEval 2021
tasks. The numbers represent normalized word counts for words classified as
containing a specific sentiment or emotion for data from the Danish dataset.

##  4 Results

For Danish, using Nordic BERT for the final submission, we obtained an
accuracy of 92.38% and F1-score of 81.18% (Figure 1 a). The confusion matrix
shows that only 9 out of 34 offensive tweets (OFF) were misclassified, and 16
out of 294 not offensive tweets (NOT) were wrongly classified as offensive.
Regarding the other languages, for Greek we obtained an F1-score of 82.6%,
77.2% for Turkish, and 73.1% for Arabic.

Figure 1: Official results of the LT@Helsinki Team.
Left: Danish, sub-task A, Right: Sub-task C (English)

In sub-task C our submission ranked second with an accuracy of 79.74% and
F1-score of 66.99%. The relatively low F1-score is due to the high numbers of
misclassifications of the minority class (only 33 OTH instances were correctly
classified out of 82) (see Figure 1 b). Despite our efforts in balancing the
dataset, it seems like the skewed class distribution inevitably added bias to
the model.

It is impressive how good BERT performed overall even though we applied very
few pre-processing steps to the data. Moreover, the results obtained on
previously unseen data indicate that none of our models were heavily
overfitted.

##  5 Conclusions

Our scores on evaluation data show that models pre-trained on general corpora
can obtain competitive results when there is very little data available, as
was the case of the Danish sub-task where we obtained the best results. BERT
performed well on both tasks, but there is still room for improvement. In the
future, we intend to experiment with those languages that we could not focus
on this time. It should be noted that multilingual BERT works best with
languages similar to English [Pires et al., 2019] so it is very likely that
the other languages would have benefited from the use of language-specific
models even more than Danish as Danish is the language most similar to English
in comparison to Greek, Turkish, and Arabic.

A more thorough comparison between multilingual and language-specific BERT
would yield more definitive answers as to whether language-specific is always
best or whether that is language-dependent. Further examination and use of
emotion-laden words in offensive texts could also help with the detection of
offensive texts. Examining other baseline methods and combining them into an
ensemble model with majority voting is another approach to consider for future
work. With regards to the class imbalance problem, other techniques such as
adjusting the classification thresholds or enriching the dataset with back and
forth machine translation of the minority class could prove to be useful.

## References

  * [Burnap and Williams, 2015]  Pete Burnap and Matthew L Williams.  2015\.  Cyber hate speech on twitter: An application of machine classification and statistical modeling for policy and decision making.  Policy & Internet, 7(2):223–242.
  * [Çöltekin, 2020]  Çağrı Çöltekin.  2020\.  A Corpus of Turkish Offensive Language on Social Media.  In Proceedings of the 12th International Conference on Language Resources and Evaluation. ELRA.
  * [Davidson et al., 2017]  Thomas Davidson, Dana Warmsley, Michael W. Macy, and Ingmar Weber.  2017\.  Automated hate speech detection and the problem of offensive language.  In ICWSM.
  * [Del Vigna et al., 2017]  Fabio Del Vigna, Andrea Cimino, Felice Dell’Orletta, Marinella Petrocchi, and Maurizio Tesconi.  2017\.  Hate me, hate me not: Hate speech detection on facebook.  In Proceedings of the First Italian Conference on Cybersecurity (ITASEC17), pages 86–95.
  * [Devlin et al., 2019]  Jacob Devlin, Ming-Wei Chang, Kenton Lee, and Kristina Toutanova.  2019\.  Bert: Pre-training of deep bidirectional transformers for language understanding.  In NAACL-HLT.
  * [Gambäck and Sikdar, 2017]  Björn Gambäck and Utpal Kumar Sikdar.  2017\.  Using convolutional neural networks to classify hate-speech.  In Proceedings of the first workshop on abusive language online, pages 85–90.
  * [Kumar et al., 2018a]  Ritesh Kumar, Atul Kr. Ojha, Shervin Malmasi, and Marcos Zampieri.  2018a.  Benchmarking aggression identification in social media.  In TRAC@COLING 2018.
  * [Kumar et al., 2018b]  Ritesh Kumar, Atul Kr Ojha, Marcos Zampieri, and Shervin Malmasi.  2018b.  Proceedings of the first workshop on trolling, aggression and cyberbullying (trac-2018).  In Proceedings of the First Workshop on Trolling, Aggression and Cyberbullying (TRAC-2018).
  * [Lefever et al., 2018]  Els Lefever, Bart Desmet, and Guy De Pauw.  2018\.  Ta-cos 2018: 2nd workshop on text analytics for cybersecurity and online safety: Proceedings.  In TA-COS 2018–2nd Workshop on Text Analytics for Cybersecurity and Online Safety, collocated with LREC 2018, 11th edition of the Language Resources and Evaluation Conference. European Language Resources Association (ELRA).
  * [Liu et al., 2019]  Ping Liu, Wen Li, and Liang Zou.  2019\.  Nuli at semeval-2019 task 6: transfer learning for offensive language detection using bidirectional transformers.  In Proceedings of the 13th International Workshop on Semantic Evaluation, pages 87–91.
  * [Malmasi and Zampieri, 2018]  Shervin Malmasi and Marcos Zampieri.  2018\.  Challenges in discriminating profanity from hate speech.  Journal of Experimental & Theoretical Artificial Intelligence, 30(2):187–202.
  * [Mathew et al., 2019]  Binny Mathew, Ritam Dutt, Pawan Goyal, and Animesh Mukherjee.  2019\.  Spread of hate speech in online social media.  In Proceedings of the 10th ACM Conference on Web Science, WebSci ’19, page 173–182, New York, NY, USA. Association for Computing Machinery.
  * [Mohammad and Turney, 2013]  Saif M. Mohammad and Peter D. Turney.  2013\.  Crowdsourcing a word-emotion association lexicon.  Computational Intelligence, 29(3):436–465.
  * [Mubarak et al., 2020]  Hamdy Mubarak, Ammar Rashed, Kareem Darwish, Younes Samih, and Ahmed Abdelali.  2020\.  Arabic offensive language on twitter: Analysis and experiments.  arXiv preprint arXiv:2004.02192.
  * [Nikolov and Radivchev, 2019]  Alex Nikolov and Victor Radivchev.  2019\.  Nikolov-radivchev at semeval-2019 task 6: Offensive tweet classification with bert and ensembles.  In Proceedings of the 13th International Workshop on Semantic Evaluation, pages 691–695.
  * [Novak et al., 2015]  Petra Kralj Novak, Jasmina Smailović, Borut Sluban, and Igor Mozetič.  2015\.  Sentiment of emojis.  PloS one, 10(12):e0144296.
  * [Oksanen et al., 2014]  Atte Oksanen, James Hawdon, Emma Holkeri, Matti Näsi, and Pekka Räsänen.  2014\.  Exposure to online hate among young social media users.  Sociological studies of children & youth, 18(1):253–273.
  * [Pedregosa et al., 2011]  Fabian Pedregosa, Gaël Varoquaux, Alexandre Gramfort, Vincent Michel, Bertrand Thirion, Olivier Grisel, Mathieu Blondel, Peter Prettenhofer, Ron Weiss, Vincent Dubourg, et al.  2011\.  Scikit-learn: Machine learning in python.  Journal of machine learning research, 12(Oct):2825–2830.
  * [Pires et al., 2019]  Telmo Pires, Eva Schlinger, and Dan Garrette.  2019\.  How multilingual is multilingual bert?  arXiv preprint arXiv:1906.01502.
  * [Pitenis et al., 2020]  Zeses Pitenis, Marcos Zampieri, and Tharindu Ranasinghe.  2020\.  Offensive Language Identification in Greek.  In Proceedings of the 12th Language Resources and Evaluation Conference. ELRA.
  * [Rosenthal et al., 2020]  Sara Rosenthal, Pepa Atanasova, Georgi Karadzhov, Marcos Zampieri, and Preslav Nakov.  2020\.  A Large-Scale Semi-Supervised Dataset for Offensive Language Identification.  In arxiv.
  * [Sigurbergsson and Derczynski, 2019]  Gudbjartur Ingi Sigurbergsson and Leon Derczynski.  2019\.  Offensive language and hate speech detection for danish.  arXiv preprint arXiv:1908.04531.
  * [Vaswani et al., 2017]  Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N Gomez, Łukasz Kaiser, and Illia Polosukhin.  2017\.  Attention is all you need.  In Advances in neural information processing systems, pages 5998–6008.
  * [Waseem et al., 2017]  Zeerak Waseem, Wendy Hui Kyong Chung, Dirk Hovy, and Joel Tetreault.  2017\.  Proceedings of the first workshop on abusive language online.  In Proceedings of the First Workshop on Abusive Language Online.
  * [Wiegand et al., 2018]  Michael Wiegand, Melanie Siegel, and Josef Ruppenhofer.  2018\.  Overview of the GermEval 2018 Shared Task on the Identification of Offensive Language.  In Proceedings of GermEval.
  * [Zampieri et al., 2019a]  Marcos Zampieri, Shervin Malmasi, Preslav Nakov, Sara Rosenthal, Noura Farra, and Ritesh Kumar.  2019a.  Predicting the type and target of offensive posts in social media.  arXiv preprint arXiv:1902.09666.
  * [Zampieri et al., 2019b]  Marcos Zampieri, Shervin Malmasi, Preslav Nakov, Sara Rosenthal, Noura Farra, and Ritesh Kumar.  2019b.  SemEval-2019 Task 6: Identifying and Categorizing Offensive Language in Social Media (OffensEval).  In Proceedings of The 13th International Workshop on Semantic Evaluation (SemEval).
  * [Zampieri et al., 2020]  Marcos Zampieri, Preslav Nakov, Sara Rosenthal, Pepa Atanasova, Georgi Karadzhov, Hamdy Mubarak, Leon Derczynski, Zeses Pitenis, and Çağrı Çöltekin.  2020\.  SemEval-2020 Task 12: Multilingual Offensive Language Identification in Social Media (OffensEval 2020).  In Proceedings of SemEval.
  * [Zhang et al., 2018]  Ziqi Zhang, David Robinson, and Jonathan Tepper.  2018\.  Detecting hate speech on twitter using a convolution-gru based deep neural network.  In European semantic web conference, pages 745–760. Springer.
  * [Öhman et al., 2016]  Emily Öhman, Timo Honkela, and Jörg Tiedemann.  2016\.  The challenges of multi-dimensional sentiment analysis across languages.  In Proceedings of the Workshop on Computational Modeling of People’s Opinions, Personality, and Emotions in Social Media (PEOPLES), pages 138–142.