text:
\section{Introduction}
\label{sec:intro}
\vspace{-0.1cm}
The Standard Model of particle physics has been observed with great
precision at the experiments. The last missing piece, the Higgs boson, was discovered
by ATLAS \cite{ATLAS:2012yve} and CMS \cite{CMS:2012qbp} at the Large Hadron Collider (LHC) in 2012. However, it has some shortcomings -- it is not possible to account for
the small neutrino masses or the existence of dark matter in the SM, for instance.
It is also difficult to explain the current measurement of muon (g-2) \cite{Muong-2:2021ojo}
and the outcomes of flavor experiments \cite{Zyla:2020zbs} within the SM framework.

The small neutrino masses can be achieved
by introducing exotic particles at a high scale via the effective
dimension-five Weinberg operator at tree level vis-a-vis the seesaw mechanism.
These extra particles correspond to
a heavy fermion singlet, a scalar triplet and a fermion triplet in type I, II and III seesaw mechanisms,
respectively \cite{Mohapatra:1979ia,Konetschny:1977bn,Magg:1980ut,
Schechter:1980gr,Cheng:1980qt,Lazarides:1980nt,Mohapatra:1980yp,Foot:1988aq}.
However, there are other models, where the exotic particle content
involves one or more larger multiplets of scalars and fermions
together. These models \cite{Liao:2010cc,Babu:2009aq,Tavartkiladze:2001by,
Kumericki:2012bh,Mohapatra:1974hk,Senjanovic:1975rk,KumarAgarwalla:2018nrn,
Ko:2015uma,Cai:2011qr,Cirelli:2009uv} not only explain the small neutrino masses but
also provide an explanation for the muon (g-2) data, flavor anomalies
while also predicting a suitable dark matter (DM)
candidate in some cases. Such models also predict exotic signatures at the colliders.

A good example of such models are cascade seesaw-like
scenarios
\cite{Liao:2010cc,Babu:2009aq,Tavartkiladze:2001by,Kumericki:2012bh}, where the neutrino mass is generated via a higher
dimensional $(5+4n)$ operator, where $n=1$
is the minimal scenario with three generations of Majorana quintuplets $\Sigma_R$, with
hypercharge $Y=0$ and a scalar quadruplet $\Phi$, with hypercharge $Y=-1$.
Another example is the left-right symmetric (LRS) framework \cite{Mohapatra:1974hk,Senjanovic:1975rk}
with an $SU(2)_R$ quintuplet where the gauge group is extended
to $SU(3)_C\times SU(2)_L\times SU(2)_R\times U(1)_{B-L}$.
In these models, the presence of a right-handed neutrino
in the particle spectrum is essentially governed by the gauge
structure and hence, naturally explains the origin of light
neutrino masses. Further, the neutral component of the quintuplet
can be a DM candidate~\cite{KumarAgarwalla:2018nrn,Ko:2015uma}.
In models such as R$\nu$MDM \cite{Kumericki:2012bh, Cai:2011qr},
purely radiative neutrino masses are generated while
also providing a viable DM candidate \cite{Cirelli:2009uv}.
Models with both exotic fermions and charged scalars can also be motivated from
the little Higgs scenarios \cite{Low:2002ws}, where the global
symmetry is $SU(6)/Sp(6)$. These models successfully explain the flavor anomalies \cite{Kumar:2020yco} and
the signatures can be studied at the colliders. Multicharged
fermions also appear in models of warped extra dimensions and
some other models as shown in Ref: \cite{Eichten:1983hw,Cabibbo:1983bk,Cirelli:2005uq,
delAguila:2008pw,Ma:2014zda,Chakrabarty:2020jro,delAguila:2010es,Wang:2017sxx}.

In this paper, we are motivated by a model with one quartet and one
septet scalar fields and quintuplet Majorana fermions (3 copies)
\cite{Nomura:2017abu}. The interactions between the SM $SU(2)_L$
lepton doublets and these large multiplets induce neutrino masses which
are suppressed by small VEV's of the quartet and/or the septet and also by
the inverse of the quintuplet fermion mass ($\sim$ TeV).
As a result, the scale of the neutrino Yukawa coupling can be reached to
less than $\mathcal{O}(1)$ \cite{Nomura:2017abu,Kumar:2019tat}. Further, this type of model is safe from any
quantum anomaly, given the zero $U(1)_Y$ charge of these quintuplet fermions.
There is no contribution to $SU(2)$ gauge anomaly as well.

The contents of the quintuplet are doubly charged fermions,
singly charged fermions and a neutral fermion\cite{Nomura:2017abu}.
Signatures of quintuplet fermions at LHC have been studied in Ref:\cite{Yu:2015pwa,Chen:2013xpa}.
These fermions are also good candidates for exotic particle search in
the future collider experiments \cite{Ozansoy:2019kdq,Guo:2016hjt}.
In most of the phenomenological studies involving the quintuplet fermions,
they decay directly into the SM particles. Even in scenarios as in
Ref: \cite{Kumericki:2012bh}, where interactions
between the quintuplet fermions and the scalar multiplets are allowed,
the fermions can not decay into the scalars as the scalars are
slightly heavier than the fermions. However, in our scenario\cite{Kumar:2019tat}, a small
mass difference between the scalar multiplets
and the fermionic quintuplet is naturally implied
from muon $(g-2)$ and in such a way that the quintuplet fermions are heavier than the scalars.

In this paper, we study the signatures of the singly and the doubly charged
quintuplet fermions ($\Sigma_R$) at the hadron collider ($pp$ collider) and linear colliders
such as $e^+e^-$ colliders,
with each having its own advantages. Although the LHC has a much higher energy reach,
the $e^+e^-$ colliders provide a cleaner environment for distinguishing the
signal from the background \cite{CLICDetector:2013tfe} and are more suitable for precision
measurements \cite{Lukic:2014mha}. With many $e^+e^-$ colliders, such as FCC-ee, ILC and CLIC in
development stage, it would be interesting to study the discovery potential
of the exotic quintuplet fermions at these colliders
\footnote{This model also offers interesting
signatures for the charged scalars which we do not address here.}.

Cases studied so far include multilepton and multijet
signatures of the quintuplet fermions, both at the LHC \cite{Yu:2015pwa, Kumericki:2012bh}
and the linear colliders \cite{Zeng:2016tmw,Ozansoy:2019kdq,Guo:2016hjt}.
As already stated, the quintuplet fermions in the scenario considered here decay into the SM particles via the scalars which leads to final states containing a large number of
leptons and jets. Such signals for quintuplet fermions have not been studied before. Even though it is difficult to reconstruct the masses of the fermion quintuplets or the scalars, given these many particle final states, we show that by carefully
choosing the final states from amongst the many possibilities, it is indeed possible to reconstruct both the masses. The masses of the exotic fermions,
such as vectorlike quarks and leptons, are constrained to be more than $\sim$ 1 TeV
\cite{Rappoccio:2018qxp, Buckley:2020wzk} and $\geq$ 740 GeV \cite{CMS:2018cgi} respectively.
However, for this model, we choose to explore masses much below 1 TeV as well as masses
larger than 1 TeV, considering the nonstandard decay modes
of the quintuplet fermions.
For the singly charged scalar, the direct search limit from LEP is 80 GeV \cite{ALEPH:2013htx}.
However, the limit on the production cross-section of the singly charged scalar
as a function of its mass is given in Ref:\cite{CMS:2019bfg}.
Multilepton states where the doubly charged scalar has decayed into two
same sign leptons are already studied in Ref:\cite{ATLAS:2017xqs,Chakrabarty:2020jro,ATLAS:2021jol}.
The lower mass limit ranges from about 230-870 GeV. However,
these studies are based on several assumptions which
do not apply to our case. In the analysis done here, the mass of the scalars is chosen to be well below
the mass of the quintuplet fermions to facilitate the decay of the quintuplet fermion via the scalars.

The rest of the paper is arranged as follows: The model is described in section \ref{sec:model}.
Collider signatures at the LHC and linear collider are described in
sections \ref{sec:pp} and \ref{sec:ee}, respectively.
We discuss the results and conclude in section \ref{sec:outlook}.
\section{Model}
\label{sec:model}
\vspace{-0.1cm}

We consider a simple scenario based on the models proposed in Ref:\cite{Nomura:2017abu,Kumar:2019tat}.
In the model of Ref:\cite{Nomura:2017abu},
the interaction among the SM $SU(2)_L$ lepton doublets and scalar multiplets, quartet
($\phi_4$, $Y=1/2$) and septet ($\phi_7$, $Y=1$),
can induce neutrino masses, while preserving the $\rho$ parameter. In order to achieve that
it is also necessary to introduce a quintuplet Majorana fermion $\Sigma_R$, with $Y=0$.
The neutrino masses are suppressed by the small VEVs of the quartet or septet and an
inverse of the quintuplet fermion mass, which explains the smallness of the neutrino mass
while also relaxing the Yukawa hierarchies.
In order to make the generation of neutrino mass more natural,
additional quintet scalar field ($\phi_5$, $Y=0$) with a
non-zero VEV can be introduced \cite{Kumar:2019tat}. As a result, tree level neutrino mass is forbidden and the
quartet scalar is an inert scalar while neutrino mass is generated at one-loop level.
The model is further constrainted from lepton flavor violation (LFV) and the muon anomalous
magnetic moment ($\Delta a_{\mu}$).

In these models, the interesting point to note is that the quintuplet fermion
$\Sigma_R$ does not decay to standard model particles directly. Instead, it happens via
its interaction with the charged scalars. Our objective is to study this particular scenario. Hence,
we choose a simplistic scenario with quintuplet fermion, $\Sigma_R$ ($Y=0$) and the
quartet scalar, $\phi_4$ ($Y=1/2$), along with their components \footnote {Even in the models with
more than one scalar multiplet, the components mix during mass diagonalization and
charged and neutral scalar mass eigenstates are obtained.}

\vspace{-0.6cm}

\begin{align}
& \Phi_4 = \left( \varphi^{++}, \varphi^{+}_2, \varphi^{0}, \varphi^{-}_1 \right)^T, \nonumber \\
& \Sigma_R = \left[ \Sigma_1^{++}, \Sigma_1^{+}, \Sigma^{0}, \Sigma_2^{+}, \Sigma_2^{++} \right]_R^T,
\end{align}

($\Sigma_1^{\pm}$, $\Sigma_1^{\pm \pm}$) and ($\Sigma_2^{\pm}$,$\Sigma_2^{\pm \pm}$)
are combined to make singly and doubly charged Dirac fermions, which we denote as $\Sigma^{\pm}$
and $\Sigma^{\pm \pm}$ respectively while $\Sigma^0_R$
remain a neutral Majorana fermion. The masses of each component are given by $M_\Sigma$ at
the tree level where mixing between the SM leptons is negligibly small.
The Yukawa interaction can be written as,

\vspace{-0.6cm}

\begin{align}
-\mathcal{L}_{Y}
&=
(y_{\ell})_{ii} \overline L_{L_i} H e_{R_i} +(y_{\nu})_{ij} [\overline L_{L_i} \tilde\Phi_4 \Sigma_{R_j} ]
 +  (M_{R})_i [\overline \Sigma^c_{R_i} \Sigma_{R_i}] + {\rm h.c.},
 \end{align}

The components of the quintuplet fermion can decay into quartet scalars and SM
leptons via the Yukawa interaction given by the second term of the above equation as,
\begin{align}
-{\cal L}_{yuk}& \supset (y_{\nu})_{ij} [\overline L_{L_i} \tilde\Phi_4 \Sigma_{R_j} ]+{\rm h.c.} \nonumber \\
&= (y_\nu)_{ij}
\left[
\overline\nu_{L_i} \left(\frac{1}{\sqrt2}\Sigma_{R_j}^0 \varphi^{0*} +\frac{\sqrt3}{2} \Sigma_{1R_j}^+ \varphi_1^- +  \frac{1}{2}\Sigma_{2R_j}^+ \varphi_{2}^- +  \Sigma_{1R_j}^{++} \varphi^{--}  \right)\right.
\nonumber\\
&\left. +
\overline\ell_{L_i} \left(\frac{1}{\sqrt2}\Sigma_{R_j}^0 \varphi_1^-  + \frac12 \Sigma_{1R_j}^+ \varphi^{--} +  \frac{\sqrt3}{2}\Sigma_{2R_j}^- \varphi^{0*} +  \Sigma_{2R_j}^{--} \varphi_{2}^{+}  \right)\right] +{\rm h.c.}
\nonumber\\
\label{eq:yukawa}
\end{align}
The coupling $y_{\nu}$, as given in Eq.~\ref{eq:yukawa}, can be constrained
from observables such as $\nu$-mass, $\Delta a_{\mu}$ and flavor observables.
Following the benchmark points in \cite{Nomura:2017abu, Kumar:2019tat}, this
coupling varies between (0.001--1), depending on the scalar particle content
of the model. As we are studying a more general scenario, for the phenomenological purpose,
we choose a very conservative limit of 0.1 for $y_{\nu}$.

The components of the quintuplet can be produced via gauge interactions given by,
\begin{align}
\overline \Sigma_R \gamma^\mu i D_\mu \Sigma_R \supset  &  \overline \Sigma^{++} \gamma^\mu \left( 2 e A_\mu + 2 g c_W  Z_\mu  \right) \Sigma^{++} +
\overline \Sigma^{+} \gamma^\mu \left( e A_\mu + g c_W  Z_\mu  \right) \Sigma^{+} \nonumber \\
&  - \sqrt{2} g   \overline \Sigma^{++} \gamma^\mu W_\mu^+  \Sigma^{+} - \sqrt{3} g   \overline \Sigma^{+} \gamma^\mu W_\mu^+ \Sigma_R^{0} - \frac{\sqrt{5} g}{\sqrt{2}}   \overline \Sigma^{+} \gamma^\mu W_\mu^+ \Sigma^{0c}_R \nonumber \\
& - \sqrt{2} g   \overline \Sigma^{+} \gamma^\mu W_\mu^-  \Sigma^{++} - \sqrt{3} g   \overline \Sigma^{0}_R \gamma^\mu W_\mu^- \Sigma^{+} - \frac{\sqrt{5} g}{\sqrt{2}}   \overline \Sigma^{0 c}_R \gamma^\mu W_\mu^- \Sigma^{+},
\label{eq:gauge_sigma}
\end{align}
where $s_W(c_W) = \sin \theta_W (\cos \theta_W)$ with the Weinberg angle $\theta_W$\footnote{The components of $\phi_4$ can also be produced at the collider via the gauge interaction,
as shown in Ref:\cite{Nomura:2017abu}.}.

The relevant gauge
interactions associated with $\phi_4$ can be obtained from the following kinetic term
\begin{equation}
|D_\mu \Phi_4|^2 \supset \sqrt{\frac{3}{2}} v_4 W^\pm W^\pm \varphi^{\mp \mp} +
\frac{g^2 v_4}{c_W} \left[ s_W^2 Z_\mu W^{+ \mu} \varphi^-_2 + \frac{\sqrt{3}}{2}(2 - s_W^2 )    Z_\mu W^{+ \mu } \varphi^-_1 + c.c. \right].
\label{eq:gauge_phi}
\end{equation}

\begin{figure}[!t]
  \centering

\begin{tikzpicture}[line width=1.4 pt, scale=1,every node/.style={scale=1.0}]
\draw[fermion,black,thin] (-4.5,1) --(-3.5,0);
\draw[fermion,black,thin] (-3.5,0) --(-4.5,-1);
\draw[vector,red,thin] (-3.5,0) --(-2,0);
\draw[fermion,black,thin] (-1,-1) --(-2,0);
\draw[fermion,black,thin] (-2,0) --(-1,1);

\node at (-4.75,1.25) {$q$};
\node at (-4.75,-1.25) {$\overline{q}$};
\node at (-2.75,0.5) {$\gamma/Z$};
\node at (-0.75,-1.25) {$\Sigma^{++}$};
\node at (-0.75,1.25) {$\Sigma^{--}$};

\draw[vector,red,thin] (1,1) --(2.75,1);
\draw[fermion,black,thick] (2.75,1) --(2.75,-1);
\draw[vector,red,thin] (2.75,-1) --(1,-1);
\draw[fermion,black,thin] (4.5,-1) --(2.75,-1);
\draw[fermion,black,thin] (2.75,1) --(4.5,1);

\node at (0.75,1) {$q$};
\node at (0.75,-1) {$\overline{q}$};
\node at (3.25,0) {$\Sigma^{++}$};
\node at (5,-1) {$\Sigma^{--}$};
\node at (5,1) {$\Sigma^{++}$};

\end{tikzpicture}
  \caption{Feynman diagrams for the production of doubly-charged fermions ($\Sigma^{\pm\pm}$) at $pp$-collider.}
\label{fig:feyn_diag_pp}
\end{figure}

Here, we have neglected the mixings among the individual components of the multiplets and we choose the components of
the quintuplet fermion as well as scalar multiplets to be degenerate, which we denote as
$M_\Sigma$ and $M_\phi$ respectively\cite{Nomura:2017abu}. We also consider
$M_\Sigma>M_\phi$, which is naturally implied form the muon anomalous magnetic moment
measurement, as shown in \cite{Kumar:2019tat}. For the complete Lagrangian
we refer to  Ref:\cite{Nomura:2017abu}.

Further, considering the interactions in Eq. \ref{eq:yukawa},
$\Sigma^{\pm}$ and $\Sigma^{\pm\pm}$ can decay via the following modes,
\begin{eqnarray}
\Sigma^{\pm} &\rightarrow &\phi_2^{\pm} \nu (\overline{\nu}) \nonumber \\
\Sigma^{\pm} &\rightarrow &\phi^{\pm \pm} l^{\mp} \nonumber \\
\Sigma^{\pm} &\rightarrow &\phi^{0} l^{\pm}
\label{eqn:decayss1}
\end{eqnarray}
\begin{eqnarray}
\Sigma^{\pm \pm} &{\rightarrow}&\phi^{\pm \pm} \nu (\overline{\nu}) \nonumber \\
\Sigma^{\pm \pm} &{\rightarrow}&\phi^{\pm} l^{\pm}
\label{eqn:decayss2}
\end{eqnarray}
The branching ratios are assumed to be same in each decay mode of $\Sigma^{\pm}$ and $\Sigma^{\pm\pm}$.
Interactions in Eq. \ref{eq:gauge_phi} allow the decays of the charged scalars into
SM gauge bosons $viz$,
\begin{eqnarray}
\phi_2^{\pm}(\phi_1^{\pm})&\rightarrow & W^{\pm} Z \nonumber \\
\phi^{\pm \pm}~~~~ &\rightarrow & W^{\pm} W^{\pm} \nonumber \\
\phi^{0}~~~~~~ &\rightarrow & W^{+} W^{-}.
\label{eqn:phidecay}
\end{eqnarray}
In the following sections, we  perform a collider study of the $\Sigma^{\pm}$ and $\Sigma^{\pm\pm}$ when they
decay via the charged scalars \footnote {We denote $\phi_2^{\pm}= \phi^{\pm}$ in the following sections.}.
Even though there are many possible production modes, \cite{Yu:2015pwa}, we study the pair
production of $\Sigma^{\pm}$ and $\Sigma^{\pm\pm}$. This is because, it will be easier to reconstruct the masses
of the quintuplets if they are produced in pairs. As can be seen from the decay modes,
the final states will have a rich collection of leptons and jets. The phenomenology of these alternative
signatures is what we study next.

\vspace{-0.1cm}

\section{Phenomenology at the $p p$ Collider}
\label{sec:pp}
\vspace{-0.1cm}

In this section, we discuss the collider physics of the quintuplet fermions at the LHC.
They can be produced in $pp$ collisions through $s$ and $t$ channel processes via
$Z/\gamma$ and $\Sigma^{\pm}$ (or $\Sigma^{\pm \pm}$), respectively.
The cross-sections for the production of the singly charged fermions,
$p p \rightarrow \Sigma^{+} \Sigma^{-}$,
are smaller in comparison to those for the doubly charged fermions,
$p p \rightarrow \Sigma^{++} \Sigma^{--}$,
as shown in Ref:\cite{Nomura:2017abu}.
Moreover, we have included the photon-photon fusion process, hence the
matrix element squared of pair productions are enhanced by a factor of $(Q)^4$,
where $Q$ is the charge of the fermion. We found the
signal to background ratio to be small for the singly charged fermions.
Hence, we choose to study the pair production of only doubly charged
fermions at the LHC and we show the Feynman diagrams in Fig:~\ref{fig:feyn_diag_pp}.
In Fig:~\ref{fig:ppssa}, we have shown the cross-section for the pair-production,
$p p \rightarrow \Sigma^{++} \Sigma^{--}$,
for different values of $\sqrt s$ at LHC where $p = q, \overline q, \gamma$. For comparison,
cross-section for pair-production of the singly charged fermion at $\sqrt s=14$ TeV
is also shown by the dotted curve in blue.

The inclusion of the photon PDF increases the signal
cross-section significantly. Moreover,
inclusion of photon PDF is important for the consistency of the calculation
as the other PDF's are determined up to NNLO in QCD. We would like to note that, in
view of the above, NNPDF \cite{NNPDF:2014otw,Ball:2013hta}, MRST \cite{Martin:2004dh}
and CTEQ \cite{Schmidt:2015zda} already include photon PDF in their definitions.
In order to compute the cross-sections and generate events at the LHC,
we incorporate the model Lagrangian in FeynRules (v2.3.13)
\cite{Alloul:2013bka,Christensen:2008py}. Using FeynRules, we generate the
model file for MadGraph5\_aMC@NLO (v2.2.1)~\cite{Alwall:2014hca}. For the
cross-sections, we use {the} NNPDF23LO1 parton distributions~\cite{Ball:2012cx}
with the factorization and renormalization scales at the
central $m_T^2$ scale after $k_T$-clustering of the event.

\begin{figure}[h!]
\begin{center}
\includegraphics[width=7cm,height=6cm]{ppssa.pdf}
\caption{Signal cross-section for the process, $p p \rightarrow \Sigma^{++} \Sigma^{--}$, where $p = q, \overline q, \gamma$
as a function of $M_\Sigma$ at different  $\sqrt s$ is shown by the solid lines. The dotted line represents the process, $p p \rightarrow \Sigma^{+} \Sigma^{-}$ at 14 TeV.}
\label{fig:ppssa}
\end{center}
\end{figure}

\vspace{-0.7cm}

\subsection{Signal}

\vspace{-0.4cm}

Once pair-produced at the LHC, the decays of the doubly charged fermions produces the following states,
\begin{equation}
\Sigma^{++} \rightarrow \phi^{++} \nu     \rightarrow (W^{+} W^{+}) \nu; ~~~
\Sigma^{--} \rightarrow \phi^{--} \overline{\nu} \rightarrow  (W^{-} W^{-}) \overline{\nu} \nonumber
\end{equation}
\begin{equation}
\Sigma^{++} \rightarrow \phi^{++} \nu     \rightarrow (W^{+} W^{+}) \nu;~~~ \\
\Sigma^{--} \rightarrow \phi^{-} l^{-}       \rightarrow  (W^{-} Z) l^{-} \nonumber
\end{equation}
\begin{equation}
\Sigma^{++} \rightarrow \phi^{+} l^{+}       \rightarrow  (W^{+} Z) l^{+};~~~ \\
\Sigma^{--} \rightarrow \phi^{-} l^{-}       \rightarrow  (W^{-} Z) l^{-}
\end{equation}
with conjugate processes included in each case.
The branching ratio in each case is assumed to be the same, as discussed in the previous section.
This gives rise to final states comprising of
a number of leptons, jets, and missing energy,
resulting in various multilepton, multijet and mixed states.
After carefully analyzing each of them on
the basis of performance over SM backgrounds and
mass reconstruction of the doubly
charged fermions, we decide upon two channels:
\begin{itemize}
\item Channel I: $\geq 4 \ell$ channel with two pairs of Same Sign (SS) leptons $(l^+l^+),(l^-l^-)$ + MET, where both the pairs are oppositely charged,
\item Channel II: $\geq 3 \ell$ + 2 jets channel with at least one pair of SS lepton $(l^\pm l^\pm)$+ One isolated lepton  $(l^\mp)$ + MET,
\end{itemize}
Here, $\ell = e, \mu, \tau$.
We also check the efficiency for the the channel $4\ell + 1/2$ jets,
but the efficiency turns out to be less compared to channels I and II.
For the $4\ell + 3/4$ jets channel, the signal and background cross-sections,
both are significantly low.
Even though it is possible to obtain a better significance
in this channel, the number of signal events to be observed at $3000 fb^{-1}$
integrated luminosity is less than 10. Hence, we do not pursue this channel.
Also, in purely leptonic channel, the leptons in the final
states are either coming from the decay of $\Sigma^{\pm\pm}$ or
from $W/Z$. Hence, the transverse mass of the leptons
can be reconstructed,
but clear mass reconstruction of the quintuplet mass is difficult
as there are different sources of MET. However, $\geq 4\ell$ channel with
the additional criteria as in
channel I, predicts a good $S/B$ ratio, and hence we study it.
The requirement of at least 4 leptons
as well as the two SS lepton pairs in channel I makes it much cleaner
compared to other channels, even though we do not put any restriction
on the number of jets. We do not impose any jet or b-jet veto
in channel-I as this will result in lesser number of signal events. Also,
we do not go beyond the 4-lepton requirement because the signal cross-section$\times$BR falls
off due to smaller branching ratio of $W$ and $Z$ into leptons.
In channel II, we include both leptons and jets in the final states
and we show that a clear mass reconstruction of $\phi^{\pm}$ or $\phi^{\pm\pm}$
\footnote{As they have the same mass.} and $\Sigma^{\pm\pm}$ is possible.

We choose five benchmark points (BP) in our study $viz$
BP1: ($M_\phi=200$ TeV, $M_\Sigma=300$ TeV)
BP2: ($M_\phi=500$ TeV, $M_\Sigma=600$ TeV), BP3:($M_\phi=700$ TeV, $M_\Sigma=800$ TeV),
BP3:($M_\phi=900$ TeV,$M_\Sigma=1000$ TeV), BP4:($M_\phi=1000$ TeV, $M_\Sigma=1200$ TeV).
We do not study the signal for $M_{\Sigma}$ larger than 1.2 TeV
due to two reasons. Firstly, the cross-section
is small at a larger mass of $\Sigma$. Secondly, when $M_\Sigma > 1.2$ TeV,
the decay products of $W$ and $Z$ bosons become colimated and the probability of
observing them as a fatjet is larger and the analysis process will
be very different. The fatjet scenario will also be ideal for 27 TeV com energy.

\vspace{-0.3cm}
\subsection{Backgrounds}

\vspace{-0.3cm}

The main backgrounds for channels I and II come
from inclusive diboson production, $VV$+jets, where $V=W,Z$.
There will also be contributions from triboson
($VVV$+jets), $HV$+jets, $t\overline t V$+jets.
The contributions from the $t\overline t VV$+ jets and
4-top backgrounds are found to be negligible.
For channel II, additionally, we get comparatively
less contribution
from $t\overline t$+ jets and $V+jets$ and
for channel I, they are negligible.
In \cite{Choudhury:2021nib} and the references therein, the cross-section of
these channels has been discussed in detail.

\vspace{-0.3cm}

\subsection{Collider Analysis}

\vspace{-0.3cm}

As a potential signature, we prefer the same sign (SS) lepton pairs
over the opposite sign (OS) leptons, due to the abundance of the former in the signal.
This is because of the decay of the quintuplets via the charged scalars, as shown earlier.
On the other hand, the SM backgrounds involving one or more than one
$Z$, are more likely to involve OS pair of leptons. Hence the signatures
involving the SS pair of leptons suffer from less SM background.
The signal and background are optimised over a set of selections, which we list
in Table \ref{tab:1}.
\begin{table}[tpb]
\begin{small}
  \begin{center}
\begin{tabular}{|c|c|c|}
\hline
   Selections & Channel I & Channel II \\
\hline
\bf S0 & $N(\ell) \geq 4$                       &  $N(\ell) \geq 3 + N(j) \geq 2$                    \\
\hline
\bf S1 & $p_T(l)> 30$ GeV                            & $p_T(l)> 30$ GeV            \\
                 & $|\eta|(l) < 2.5$                &  $|\eta|(l) < 2.5$            \\
                 & $\Delta R_{\ell \ell/j} > 0.3$    & $\Delta R_{\ell \ell/j}> 0.3$            \\
                 &                                   & $p_T(j)> 30$ GeV            \\
                 &                                  &  $|\eta|(j) < 2.5$          \\
                 &                            &  $\Delta R_{jj} > 0.3$       \\
\hline
\bf S2     &  $S_T(l) > 400 \ensuremath{\mathrm{Ge\kern -0.1em V}}$ & $S_T(l) > 200$ GeV \\
            & $M (\ell^+, \ell^+), M(\ell^-, \ell^-) > 100\ensuremath{\mathrm{Ge\kern -0.1em V}}$ &  $M (\ell^\pm, \ell^\pm) > 100\ensuremath{\mathrm{Ge\kern -0.1em V}}$        \\

            &$\Delta R(\ell_0,\ell_1) > 1.5$   &$\Delta R(\ell,\ell) > 1.5$\\
             &                &    $H_T(j) > 300$ GeV                       \\
\hline
 \bf S3     &  &  $60 \ensuremath{\mathrm{Ge\kern -0.1em V}} <M(j, j) < 120\ensuremath{\mathrm{Ge\kern -0.1em V}}$ \\
                  &                  &   $(M_\phi-100)<M(j,j,\ell)<(M_\phi+100)$ \\
\hline

\end{tabular}
\caption{\label{tab:1} Selections $S1$, $S2$ and $S3$ for channels I and II.}
\end{center}
\end{small}
\end{table}

\vspace{-0.2cm}
\paragraph{(I) \underline {$\geq 4\ell$ channel with $(l^+l^+)$ and $(l^-l^-)$ pair+MET:}}
In the multilepton searches performed for exotic particles
(vector like leptons, charged scalars etc.),
only one OS or SS pairs are identified in most of the cases, the only exception being
the searches for multicharged scalars \cite{ATLAS:2017xqs}.
This is largely due to the
standard decay modes of the BSM particles. Whereas, in channel I, we have the
requirement for two SS lepton pair, with pairwise opposite charge. In channel I,
the quintuplet fermions
decay to $Z$ and/or $W$ via the charged scalars, and most of the final state leptons
come from the decay of the $Z/W$ bosons. Hence, we do not apply any $Z/W$ veto.
We arrange the leptons ($\ell_i$) in the descending order of $pT$.
In this channel, the SS leptons appear directly
from the direct decay of $\Sigma^{\pm\pm}$ and from the decay of
$W/Z$ decay in the same decay chain.
We plot the lepton's $pT$ and the invariant mass distribution of
the SS lepton pairs for 2 BP's (BP2 and BP5) in Fig:~\ref{fig:dist1}.
The solid line represents the invariant mass distribution if at least one SS
pair is present, and the dotted line represents the same for the
additional SS pair, where these two pairs have
opposite charge. Based on the two body mass distribution, we have imposed
the selection $M (\ell^\pm, \ell^\pm) > 100\ensuremath{\mathrm{Ge\kern -0.1em V}}$.
\begin{figure}
\begin{center}
\includegraphics[width=5.5cm,height=4.5cm]{ptl600-eps-converted-to.pdf}
\includegraphics[width=5.5cm,height=4.5cm]{ptl1200-eps-converted-to.pdf}\\
\includegraphics[width=5.5cm,height=4.5cm]{mll600-eps-converted-to.pdf}
\includegraphics[width=5.5cm,height=4.5cm]{mll1200-eps-converted-to.pdf}
\caption{(Top) Distributions of $pT(l)$ and (Bottom) same sign lepton pair invariant mass distributions for $M_{\Sigma}=$ 600 GeV and 1200 GeV (BP2 and BP5) respectively in channel I at 14 TeV LHC. The solid and dashed lines correspond to the first and second pair of SS lepton pairs.}
\label{fig:dist1}
\end{center}
\end{figure}

In Fig:~\ref{fig:dist2},
we plot the sum of lepton $pT$ ($S_T (\ell)$) and
missing transverse energy (MET) distributions for the signal and the total background.
Note that, a substantial amount of
MET is present in the signal, as well as in the SM background.
Hence, we restrain from putting any cut on MET in order to get most of the signal events.
Even though the signal has a very high $S_T (\ell)$ compared to the background,
the set of selections optimise for $S_T (\ell)> 400$ GeV for the whole
signal region under consideration. For example, if we focus on the region $M_\Sigma > 1$ TeV only,
$S_T  (\ell)> 600$ GeV gives a much better $S/B$ ratio. But wwe choose to use only one value for the
$S_T  (\ell)$ selection for our BP points.
Based on the plots of the kinametic variables, we optimize the selections at the given
values in Table \ref{tab:1}. We found that
the cut on $S_T(\ell)$
is sufficient to suppress the background in channel I.
Moreover, the leptons with heighest $pT$ will have a large separation
compared to the other leptons, as they are form the separate
decay chains of the quintuplet in most of the cases.
Thus we impose a selection on these leptons by requiring
$\Delta R (\ell_0$,$\ell_1)>1.5$.
\begin{figure}
\begin{center}
\includegraphics[width=5.5cm,height=4.5cm]{met-eps-converted-to.pdf}
\includegraphics[width=5.5cm,height=4.5cm]{htlept-eps-converted-to.pdf}
\caption{(left) Transverse missing energy (MET) and (right) sum of lepton $pT$ ($S_T(\ell)$)
distribution for $M_{\Sigma}=$ 600 GeV and 1200 GeV, and
total background (shadowed region) in channel I.}
\label{fig:dist2}
\end{center}
\end{figure}
In Table~\ref{tab:2}, we summarize the effect of the
selections in channel I with $(l^+l^+)$ and $(l^-l^-)$ pair + MET.
As there is no jet veto imposed, the majority of the backgrounds come from
the diboson + jets events. Initially, this background cross section is comparably very high but after we
impose the selection $S_2$, the background reduces further.
\begin{table}[!h]
\parbox{.5\linewidth}{

~~~~~~~~~~~~~\begin{tabular}{|c|c|c|c|}
\hline
\hline
$M_{\Sigma}$ & $S1$ (fb) & $S2$ (fb) \\
\hline
BP1, 300 GeV& 2.335		&1.093	\\
BP2, 600 GeV&0.598		&0.219	\\
BP3, 800 GeV&0.196		&0.093\\
BP4, 1000 GeV&0.063		&0.032\\
BP5, 1200 GeV&0.024		&0.013\\
\hline
\hline
\end{tabular}
}
\parbox{.5\linewidth}{
\begin{tabular}{|c|c|c|c|}
\hline
\hline
Major Backgrounds &$S1$ (fb)& $S2$(fb)\\
\hline
Di-Boson+jets& 20.42 & 0.55\\
$t\overline{t}V$       &0.25&0.07 \\
Triboson &0.082 & 0.022\\
$HV$+jets &0.048& 0.021 \\
\hline
Total   &20.80&0.66 \\
\hline
\hline
\end{tabular}
}
    \caption{\label{tab:2}\em {\em (Left)} The variation of the cross-section (fb) for each of the
    BPs, as the selections are imposed at 14 TeV LHC in channel I. {\em
    (Right)} The same is shown for the background.}
\end{table}
\paragraph{(II) \underline {$\geq 3\ell$ Channel with $(l^\pm l^\pm)$ pair + $l^\mp$ + $\geq 2$ jets channel}:}
In this channel, we require the presence of at
least three leptons as well as two or more jets.
In Fig:~\ref{fig:dist3}, we show the $p_T$ distributions
of the jets, and also the sum of $pT$
for the jets ($H_T(j)$). The $pT$ distribution of the
leptons are mostly the same as channel I but
as the number of leptons in channel II
is less than channel I, the
$S_T(\ell)$ distribution peaks at a lower value
compared to channel I. In channel II,
we identify at least one SS lepton pair in a manner as stated in channel I.
The selections in channel II are summerized in Table \ref{tab:1}.
Additionally,
we find that a cut on the minimum value of $H_T(j)$ is useful
to minimise the background.

\begin{figure}
\begin{center}
\includegraphics[width=5.5cm,height=4.5cm]{ptj600-eps-converted-to.pdf}
\includegraphics[width=5.5cm,height=4.5cm]{ptj1200-eps-converted-to.pdf}\\
\includegraphics[width=5.5cm,height=4.5cm]{htlept1-eps-converted-to.pdf}
\includegraphics[width=5.5cm,height=4.5cm]{htjets-eps-converted-to.pdf}
\caption{
(Top) Distributions of jet $pT$ and (Bottom) sum of lepton $pT$ ($S_T(\ell)$) (left) and sum of jet $pT$ ($H_T(j)$) (right), for $M_{\Sigma}=$ 600 GeV and 1200 GeV, respectively, for channel II at 14 TeV LHC. The shadowed region corresponds to the total background.}
\label{fig:dist3}
\end{center}
\end{figure}

The main objectives in channel II are to
construct the three and four body
invariant mass distribution, $M(\ell \ell j)$ and
$M(\ell\ell jj)$, for the reconstruction of
$M_\phi$ and $M_\Sigma$, respectively. Note that, this is only possible
when we consider the decay of the quintuplet via the singly charged scalar.
Even though it is theoretically possible to reconstruct the mass of the quintuplet from
$M(\ell jjj)$ also, it is harder to select the exact jets for the distribution.
Hence, we consider
the case when $W$ and $Z$ decay through leptonic mode and hadronic mode respectively.
At first, we select two SS leptons in such a way that they must come from the same decay chain.
One lepton is coming from the decay of the quintuplet and the another is from the $W$. We demand
$\Delta R(\ell\ell) >1.5$ for these two leptons.
Then we select two jets coming from the the decay of the $Z$ boson,
by requiring $60 < M(j j) < 120$ GeV.
The three body mass distribution $M(\ell \ell j)$
and the four body mass distribution
$M(\ell\ell jj)$ reconstruct the masses of $\phi$ and $\Sigma^{\pm}$, which is shown
in Fig:~\ref{fig:dist4}.
\begin{figure}
\begin{center}
\includegraphics[width=5.5cm,height=4.5cm]{mljj2-eps-converted-to.pdf}
\includegraphics[width=5.5cm,height=4.5cm]{mlljj2-eps-converted-to.pdf}
\caption{(left) Three body
invariant mass $M(\ell\ell jj)$ and (right) four body invariant mass $M(\ell\ell jj)$ for $M_{\phi}=$ 500, 700, 900 GeV and $M_{\Sigma}=$ 600, 800, 1000 GeV, respectively, for channel II at 14 TeV LHC. The shadowed region corresponds to the total background.{\bf change plot}}
\label{fig:dist4}
\end{center}
\end{figure}
We select the final events with $S_3$, where
the events are required to satisfy the three body
invariant mass in the window of $M_\phi\pm 100$ GeV.
The signal and background cross-section after the cuts
are shown in Tab~\ref{tab:3}.
Clearly, the selection after $S_3$ gives a better signal
to background ratio.
\begin{table}[!h]
  \begin{center}
\begin{tabular}{|c|c|c|c|}
\hline
\hline
$M_{\Sigma}$&$S1$ & $S2$ & $S3$ \\
\hline
BP1, 300 GeV& 11.55 & 4.4 & 0.112\\
BP2, 600 GeV& 2.85 &1.025 &0.028	\\
BP3, 800 GeV& 0.84 &0.39	 &0.007\\
BP4, 1000 GeV& 0.25 &0.14	 &0.0017	\\
BP5, 1200 GeV& 0.09 &0.054 &0.0005	\\
\hline
\hline
\end{tabular}
\vskip 10pt
\begin{tabular}{|c|c|c|ccccc|}
\hline
\hline
Major Backgrounds &$S1$ & $S2$ & $S3$ (BP1)& $S3$ (BP2) & $S3$ (BP3) &$S3$ (BP4)&$S3$ (BP5)\\
\hline
Di-Boson+jets   & 45.05 &14.00 & 0.085  &5.5$\times 10^{-3}$  & 3$\times 10^{-4}$ &$6\times 10^{-5}$& $1\times 10^{-5}$ \\
$t\overline{t}V$& 10.42 &0.53  & 0.056 &1$\times 10^{-3}$    & 2$\times 10^{-4}$ &$2\times 10^{-5}$& $<10^{-5}$\\
Triboson   & 0.336 &0.013 & 0.004 & $<10^{-3}$          & $<10^{-4}$        &$<10^{-5}$       &$<10^{-6}$\\
$HV$+jets  & 1.2   &0.012 & 0.003 &$<10^{-3}$           & $<10^{-4}$        &$<10^{-5}$       &$<10^{-6}$\\
\hline
Total   & 61.06 & 14.57 &0.148& 0.0065 &0.0005 & $7\times 10^{-5}$ & $1\times 10^{-5}$\\
\hline
\hline
\end{tabular}
\end{center}
    \caption{\label{tab:3}\em {\em (Top)} The variation of the signal cross-section (fb) for each of the
    BPs, as the selections are imposed at 14 TeV LHC. {\em
    (Bottom)} The same for the various backgrounds.}
\end{table}
\paragraph{\underline {Result}}
The significance for the discovery can be described as (see Ref: \cite{Cowan:2010js,Li:1983fv,Cousins:2007yta}),

\vspace{-0.7cm}

\begin{center}

\begin{equation}
\Zdis = \left [ 2\left ((s+b) \ln \left [\frac{(s+b)(b+\sigmab^2)}{b^2+(s+b)\sigmab^2}\right ]-
         \frac {b^2}{\sigmab^2} \ln\left [1+ \frac {\sigmab^2 s}{b(b+\sigmab^2)}\right ]\right)\right]^{1/2}
         \label{eq:zdis1}
\end{equation}

\end{center}
Where $s$ and $b$ are number of signal and background events respectively, and $\Delta_b$ is the
uncertainty in the measurement of the background. If $\Delta_b = 0$,
\begin{center}
$\Zdis=\sqrt{2[(s+b)\ln(1+s/b)-s]}$
\end{center}
If $b$ is large,
\begin{center}
$\Zdis = s/\sqrt{b}$
\end{center}

Thus, if $b$ is small, $s/\sqrt{b}$ overestimates the significance. We
use $\Zdis > 5$ which corresponding to $p< 2.86 \times 10^{-7}$ for different values of
$\Delta_b$. Similarly, the significance
for exclusion is,

\vspace{-0.5cm}

\begin{eqnarray}
\Zexc&=\left [2 \left \{ s-b \ln \left (\frac{b+s+x}{2b} \right ) - \frac{b^2}{\Delta_b^2} \ln \left (\frac{b-s+x}{2b} \right ) \right \} -
(b + s - x) (1 + b/\Delta_b^2) \right ]^{1/2} \nonumber \\
&x = \sqrt{(s+b)^2 - 4 s b \Delta_b^2/(b + \Delta_b^2)}
         \label{eq:zdis2}
\end{eqnarray}

If $\Delta_b = 0$,
\begin{center}
$\Zexc = \sqrt{2(s - b \ln(1 + s/b))}$
\end{center}
For 95\% confidence level (CL) exclusion
($p = 0.05$), we use $\Zexc > 1.645$ for different
$\Delta_b $.

We calculate the significance using the formula in Eq.\ref{eq:zdis1},
 in order to account
for the uncertainty in the background, as the background is small in both the channels.
The integrated luminosity for
discovery and exclusion as a function
of the mass of the doubly charged fermion($M_\Sigma$) is shown in Fig:~\ref{fig:result1}.
The prediction in
channel I is sensitive to the uncertainty in the background,
which we have considered to be
$\sigma_B= 0$, $0.25\times b$, $0.5\times b$. Channel II is
not sensitive to $\sigma_B$ as the signal
and background cross-sections, both, are small, as given in Table~\ref{tab:3}.
We have found that channel I and II
have a good discovery potential for masses upto
850 GeV and 1025 GeV respectively at 3000 fb$^{-1}$ luminosity with $\sigma_B$=0.
In channel I, more than $3000$ fb$^{-1}$ luminosity
is required for discovery of $M_\Sigma >850$ GeV with nonzero $\sigma_B$.

Masses upto 1.05 TeV and 1.2 TeV can
be excluded with 95\% CL (corresponds to Z value =1.645)
at 3000 fb$^{-1}$ luminosity, with no background uncertainty.
In channel I, with integrates luminosity of $3000$ fb$^{-1}$,
the exclusion limit is 1 TeV and 920 GeV,
for $\sigma_B$=0.25 and 0.5 respectively.
Hence, we find that channel I and II have a good prospect
for both exclusion and discovery of the doubly
charged fermions in HL-LHC, with $3000$ fb$^{-1}$, with the added advantage of mass reconstruction
for the doubly charged fermion and the charged scalar in channel II.
\begin{figure}[tb]
\begin{center}
\includegraphics[width=6cm,height=5cm]{lumidis-eps-converted-to.pdf}
\includegraphics[width=6cm,height=5cm]{lumiexc-eps-converted-to.pdf}
\caption{The integrated luminosity for
discovery (left) and exclusion (right) as a function
of the mass of the doubly charged fermion at 14 TeV LHC. The solid, dashed and dotted lines correspond
to 0\%, 25\% and 50\% uncertainty in the total background, respectively.}
\label{fig:result1}
\end{center}
\end{figure}
\section{Phenomenology at the $e^{+}e^{-}$ Collider}
\label{sec:ee}
\vspace{-0.1cm}

We have shown in Sec:~\ref{sec:pp} that the pair production
cross-sections of the singly charged fermions ($\Sigma^{\pm}$)
are smaller compared to the doubly charged fermions ($\Sigma^{\pm\pm}$)
at p-p collision, where both are components of a fermionic quintuplet. The small cross-section makes it difficult
to observe singly charged fermions at 14 TeV LHC when
we look for the alternative signatures in our
model. Even increasing the center of mass
energy further upto 27 TeV does not solve the issue.
These singly and doubly charged fermions can also be produced
in linear colliders, such as the $e^+e^-$ collider, which in turn generate
multiple leptons and jets in the final state. Even though it is possible to
observe alternative signatures for both singly and doubly charged fermions
at  the $e^+e^-$ colliders, we restrict ourself to the case
of the singly charged fermion. The production of the doubly charged fermions
lead to more leptons and jets in the final state
than the singly charged fermions. Here, we choose to study the final states once
the singly charged fermions are produced in pair at $e^+e^-$ collider. The analysis
for the doubly charged fermions will be similar to this.
Moreover, at LHC, being a $pp$ collider, the multijet signals are complicated to study due to the heavy
QCD backgrounds. But the $e^+e^-$ collider offers a much clean environment. Hence,
the SM background for the signal involving multiple jets are
remarkably small compared to $pp$ collider.

\vspace{-0.3cm}
\subsection{Signal}
\vspace{-0.3cm}
The singly charged fermions($\Sigma^{\pm}$),
can be produced in pairs at the $e^+e^-$-collider via the
gauge couplings as described in the previous section.
The Feynman diagrams for the pair production are shown in
Fig:~\ref{fig:feyn_diag}. In general, the process proceeds through $s$-channel
via $\gamma$ and $Z$ boson exchange. But, in this particular model,
there is an extra contribution
coming from the $t$-channel diagram via the doubly charged scalar.
The cross-section due to the t-channel diagram is large compared to the
other diagram.
However, the contribution in the total cross-section is not so large
due to destructive interference between the $s$- ans the $t$-channel diagrams.
\begin{figure}[!t]
  \centering

\begin{tikzpicture}[line width=1.4 pt, scale=1,every node/.style={scale=1.0}]

\draw[fermion,black,thin] (-4.5,1) --(-3.5,0);
\draw[fermion,black,thin] (-3.5,0) --(-4.5,-1);
\draw[vector,red,thin] (-3.5,0) --(-2,0);
\draw[fermion,black,thin] (-1,-1) --(-2,0);
\draw[fermion,black,thin] (-2,0) --(-1,1);

\node at (-4.75,1.25) {$e^-$};
\node at (-4.75,-1.25) {$e^+$};
\node at (-2.75,0.5) {$\gamma/Z$};
\node at (-0.75,-1.25) {$\Sigma^+$};
\node at (-0.75,1.25) {$\Sigma^-$};

\draw[fermion,black,thin] (1,1) --(2.75,1);
\draw[scalar,blue,thick] (2.75,1) --(2.75,-1);
\draw[fermion,black,thin] (2.75,-1) --(1,-1);
\draw[fermion,black,thin] (4.5,-1) --(2.75,-1);
\draw[fermion,black,thin] (2.75,1) --(4.5,1);

\node at (0.75,1) {$e^-$};
\node at (0.75,-1) {$e^+$};
\node at (3.25,0) {$\phi^{++}$};
\node at (5,-1) {$\Sigma^+$};
\node at (5,1) {$\Sigma^-$};

\end{tikzpicture}
  \caption{Feynman diagrams for the production of singly charged quintuplet fermion at the $e^+e^-$ collider.}
\label{fig:feyn_diag}
\end{figure}
The effect of the polarization of the electron and
positron beam has been discussed in detail
in Ref: \cite{Shang:2021mgn},
and we have followed the exact same polarization of the $e^+$ and $e^-$ beam
which leads to maximum -60\% left right asymmetry ($A_{LR}$).

The production cross-sections are computed in MadGraph5\_aMC@NLO (v2.6.5) with the
normalisation and factorisation scales set at $m_Z$ and shown in Fig:~\ref{fig:prod1}.
For further study, we choose
the following benchmark points: $M_\Sigma=200$ GeV at $\sqrt s$ = 500 GeV,
$M_\Sigma=300,400$ GeV at $\sqrt s$ = 1000 GeV, and
$M_\Sigma=500, 600, 700$ GeV for $\sqrt s$ = 1500 GeV
\footnote {We did not go to masses beyond 700 GeV because it will require a
3 TeV linear collider. At large energies, the analysis will require a detailed
study of fatjets ($W/Z$) \cite{Das:2020gnt,Das:2020uer} which emerge as decay products of the
charged fermions.}.

The decays of the singly charged fermions lead to the following final states involving W/Z bosons,
\begin{figure}[h!]
\begin{center}
\includegraphics[width=5.3cm,height=5cm]{eess1_1.pdf}
\includegraphics[width=5.3cm,height=5cm]{eess1_2.pdf}
\includegraphics[width=5.3cm,height=5cm]{eess1_3.pdf}
\caption{Pair production cross-sections for the singly charged fermions of different masses
as a function of center of mass energy, at $e^+e^-$ collider.}
\label{fig:prod1}
\end{center}
\end{figure}

governed by the equations in Sec:~\ref{sec:model}. Among all final states,
the leptonic final states or final states of
leptons+jets suffer from lower effective cross-section due to small
branching ratio of $W/Z$ into leptons.
The signals with multiple jets
have the advantage over multilepton states as the branching ratio
of $W/Z$ is more into jets than leptons.
The final states involving multiple jets can have a maximum of 6 jets,
coming from the decays of $W/Z$. Here, we show a detailed
analysis of two final states:
\vspace{-0.3cm}
\begin{itemize}
\item Channel(A): One lepton $(\ell^\pm)$ + 4 jets.
\item Channel(B): Two opposite sign lepton pair $(\ell^+ \ell^-)$ + 4 jets.
\end{itemize}
\vspace{-0.3cm}
These type of signals
in $e^+e^-$ collider have a great chance for discovery due to smaller background,
which is also shown in \cite{Azuelos:1993qu}.

\vspace{-0.3cm}

\subsection{Backgrounds}

\vspace{-0.3cm}

The major backgrounds for the channels under study
get contribution from di-boson ($WW$, $ZZ$), $t\overline t$,
$t\overline V$, Triboson ($VVV$=$ZZZ$, $ZWW$) and $HZ$ production.
The variation of these
major backgrounds with $\sqrt{s}$ is already shown
in \cite{Das:2020gnt}. Along with multileptons, as
the channels under investigation
include multiple jets, we demand {\it inclusive} cross-section
of these backgrounds by producing at least two
jets in association, such as, di-boson+ 2 jets production($VVjj$),
$t\overline t$+ 2 jets and $HZ$+2 jets. The contribution from the
$\ell\ell$+ 2 jets, 4-jets and 4-top production are
found to be small. We include these backgrounds in ``others'' category.
Among all the backgrounds, the cross-section of $ZZjj$ is found to be larger.

\vspace{-0.3cm}

\subsection{Collider Analysis}

\vspace{-0.3cm}

In order to generate events,
we use MadGraph5\_aMC@NLO (v2.2.1)~\cite{Alwall:2014hca}, where the showering
and hadronization are done in a similar way as mentioned before
in the LHC part. In FastJet, the jets are reconstructed with distance parameter
$R=0.4$ using anti-$K_t$ algorithm. In Delphes, we use the Delphes ILD card \cite{Behnke:2013lya}
for detector simulation.
The signal and background events are required to pass
through selections on different kinametic distributions, as given
in Table~\ref{tab:4}. At first, we select events with basic cuts, $A1$.
Later, while selecting the single lepton or the oppositely charged lepton
pair, we make sure that it is well isolated from the jets coming from the decays
of $W/Z$ by requiring a moderate isolation cut in $A2$. We have also imposed a cut
on $M(\ell^+, \ell^-)$ in channel (B) to reduce the background further.
\begin{table}[tpb]
\begin{small}
  \begin{center}
\begin{tabular}{|c|c|c|}
\hline
   Selections & Channel (A) & Channel (B) \\
\hline
\bf A0 & $N(\ell) \geq 1$+ $N(j)\geq 4$                     &  $N(\ell) \geq 2$ + $N(j)\geq 4$                    \\
\hline
\bf A1 & $p_T(l)> 10$ GeV                            & $p_T(l)> 10$ GeV            \\
                 & $|\eta|(l) < 2.5$                &  $|\eta|(l) < 2.5$            \\
                 & $\Delta R(\ell, \ell/j) > 0.4$    & $\Delta R(\ell, \ell)> 0.4$            \\
                 &  $p_T(j)> 20$ GeV                                  & $p_T(j)> 20$ GeV            \\
                 &    $|\eta|(j) < 5.0$                                 &  $|\eta|(j) < 5.0$          \\
                 &   $\Delta R_{jj} > 0.4$                         &  $\Delta R_{jj} > 0.4$       \\
\hline
\bf A2           &     $\Delta R(\ell,j) > 1.5$       & $\Delta R(\ell,j) > 1.5$ \\
                 &      --                                                  &  $M (\ell^+, \ell^-) > 100\ensuremath{\mathrm{Ge\kern -0.1em V}}$\\
\hline

\end{tabular}
\end{center}
\caption{\label{tab:4}Selections A1 and A2 for channel (A) and channel (B).}
\end{small}
\end{table}

The signal and background cross-sections after the cuts are shown in
Table~\ref{tab:5} and Table~\ref{tab:6}, respectively.
We found the background to be small enough to give a very good signal to background ratio($S/B$),
after the initial cuts $A1$ for channel (A). For channel (B), in order to improve the $S/B$ ratio,
we have imposed further cuts in $A2$ on selected opposite sign lepton pair ($\ell^+\ell^-$), as shown in Table~\ref{tab:4}.
The requirement of exactly two leptons with opposite sign in channel (B) makes the cross
section smaller than channel (A).
The largest background
contribution comes from $t\overline{t} + jets$ due to the large cross-section.
We further check that the additional cuts on kinematic variables such as $H_T$, $S_T$ or MET
would reduce the signal efficiency effectively, hence we did not impose them.
\begin{table}[htb!]
\centering
\begin{tabular}{|c|c|c|c|c|}
\hline
\hline
\bf $\sqrt{s}$ & \bf $M_{\Sigma}$ (GeV) & \bf $\sigma$ (pb) & \bf $\sigma^{A}_{A2}$ (fb) & \bf $\sigma^{B}_{A2}$ (fb)\\
\hline
500 GeV & 200 & 0.706 & 4.45 & 0.049\\
\hline
& 200 & 0.218 & 15.70 & 0.131\\
1 TeV & 300 & 0.209  & 14.63 & 0.125\\
& 400 & 0.175 & 13.50 & 0.122\\
\hline
 & 500 & 0.089 & 7.56 & 0.107\\
1.5 TeV & 600 & 0.077 & 7.24 & 0.1001\\
 & 700 & 0.05 & 4.95 & 0.09\\
\hline
\hline
\end{tabular}
\caption{Cross-sections for the signal $e^+ e^- \rightarrow \Sigma^+ \Sigma^-$ before and after the selections.
Channel (A) corresponds to $l^{\pm}$+ 4 jets and Channel (B) corresponds to $l^{\pm}l^{\pm}$+ 4 jets.
}
\label{tab:5}
\end{table}
\begin{table}[htb!]
\centering
\begin{tabular}{|c| rrr | rrr|}
\hline
\hline
~ & \multicolumn{3}{c}{$\sigma^{A}_{A2}$ (fb)}  & \multicolumn{3}{c}{$\sigma^{B}_{A2}$ (fb)} \\
\hline
\bf Background & $\sqrt{s}=500$ GeV & 1 TeV & 1.5 TeV &  $\sqrt{s}=500$ GeV & 1 TeV & 1.5 TeV \\
\hline
$diboson+jets$ & 8.08 & 3.04 & 1.63 & 0.0 & 0.0 & 0.0\\
$t\overline{t} + jets$ & 82.5 & 24.75 & 11.25 & 1.1 & 0.33 & 0.15\\
$t\overline{t}V$ & 1.121 & 1.79 & 1.083 & 0.039 & 0.0614 & 0.037\\
$VVV$ & 2.85 & 5.0 & 3.67 & 0.0024 & 0.0035 & 0.0031\\
$HV + jets$ & 2.85 & 0.65 & 0.3 & 0.045 & 0.0 & 0.0\\
$others$ & -- & -- & -- & 0.045 & 0.0425 & 0.035\\
\hline
\hline
$Total$ & 97.4 & 34.35 & 17.94 & 1.186 & 0.437 & 0.225\\
\hline
\end{tabular}
\caption{Cross-sections for various backgrounds corresponding
to channels A ($l^{\pm}+ 4 jets$) and B ($l^{\pm}l^{\pm}+ 4 jets$) after the selections.}
\label{tab:6}
\end{table}
Even though the signal cross section in channel (B) is less, we find that it is an excellent channel to
reconstruct the invariant mass of $M_\Sigma$ from the $\ell^+ \ell^- j j$ distribution.
We show the distribution for two cases, $M_{\Sigma}$=400 and 600 GeV at $\sqrt{s}=$1 and 1.5 TeV
respectively in Fig:\ref{fig:dist}.
\begin{figure}[tb]
\begin{center}
\includegraphics[width=6.0cm,height=5.7cm]{dist1-eps-converted-to.pdf}
\includegraphics[width=6.0cm,height=5.7cm]{dist2-eps-converted-to.pdf}
\caption{Three body invariant mass
$M(\ell j j)$ for $M_{\Sigma}=$ 400 GeV (left) and $M_{\Sigma}=$ 600 GeV (right) for channel (B) at 1 TeV and 1.5 TeV $e^+e^-$ collider
respectively. The shadowed region represents the total background.}
\label{fig:dist}
\end{center}
\end{figure}

\vspace{-0.3cm}
\paragraph{\underline {Result}}
For the study in $e^+e^-$ collider, the background is larger
compared to the LHC scenario. Hence, Eq:\ref{eq:zdis1} reduces
to a simple form of $S/\sqrt{B}$. The integrated luminosity for
discovery as a function
of the mass of the singly charged fermion is shown in
Fig:\ref{fig:result2}. We find that,
the discovery potential of channel (A) is much better than channel (B), i.e,
the required integrated luminosity is less in channel (A),
for $M_\Sigma\leq 450$ GeV and $M_\Sigma\leq 750$ GeV,
at $\sqrt{s}=1$TeV and 1.5 TeV respectively.
With $\leq 20$ fb$^{-1}$ luminosity, it is possible to discover
in the region of $M_\Sigma \leq 700$ GeV with $5\sigma$ significance
in channel (A). Whereas, the required luminosity for $5\sigma$
discovery is $\leq 1000$ fb$^{-1}$, for the same in channel (B).
For 95\% exclusion limit, the entire mass region can be probed with luminosity less than
100 fb$^{-1}$ in channel (B). The required luminosity in channel (A) is very small for the same,
hence we do not plot them.
\begin{figure}[tb]
\begin{center}
\includegraphics[width=6.0cm,height=5cm]{lumi1-eps-converted-to.pdf}
\includegraphics[width=6.0cm,height=5cm]{lumi2-eps-converted-to.pdf}
\caption{The integrated luminosity for
discovery (Solid line) as a function
of the mass of the singly charged fermion ($M_\Sigma$) is shown
for channel (A) and (B) at 1 TeV (left) and 1.5 TeV (right) $e^+e^-$ collider.
The integrated luminosity for exclusion (dotted line) is also shown for channel B.}
\label{fig:result2}
\end{center}
\end{figure}

\section{Conclusion and Outlook}
\label{sec:outlook}
We have discussed the discovery potential of the singly and the doubly
charged fermions, which are components of a quintuplet, at the LHC and future
$e^+e^-$ colliders. Such a specific model, as we have considered, with quintuplet fermions
and a scalar multiplet, predicts certain signatures which require alternate search
strategies.

In the study of signatures at the LHC, we have discussed the possible multilepton
and multi(lepton+jet) signatures of the doubly charged fermions, as they have larger
cross-sections compared to that of the singly charged fermions.
For the doubly charged quintuplet fermion ($\Sigma^{\pm\pm}$),
5$\sigma$ discovery might be possible at integrated luminosity of
3000 fb$^{-1}$ at 14 TeV LHC if $M_\Sigma\leq 980$ GeV. The exclusion
limit can be extended upto 1.2 TeV with the same parameters.

On the other hand, linear colliders, such as the $e^+e^-$ collider, offer
a much cleaner environment to study the signatures associated with multiple jets. Thus,
the signals have the advantage of a larger cross-section $\times$ BR, where the $W/Z$ bosons decay
into jets. We find that the singly charged fermion ($\Sigma^\pm$) shows a great discovery potential
at the $e^+e^-$ collider, unlike the case of the LHC.
There might be a possibility of 5 $\sigma$ discovery with 1000 fb$^{-1}$
luminosity at $e^+e^-$ collider for $M_\Sigma\leq 700$ GeV.
Similar kind of final states also exist for the doubly charged quintuplet fermion
but with more leptons and jets, making the analysis much more complicated.
Thus, we will address it somewhere else.
The cross-sections for the pair production of the doubly charged fermions
at the $e^+e^-$ collider are shown in
Fig:~\ref{fig:eess2} (left and middle).
\begin{figure}[tb]
\begin{center}
\includegraphics[width=5.3cm,height=5cm]{eess2_2.pdf}
\includegraphics[width=5.3cm,height=5cm]{eess2_3.pdf}
\includegraphics[width=5.3cm,height=5cm]{aass1.pdf}
\caption{(Left and Middle) The doubly-charged quintuplet fermion production cross-section as a function of the center of mass energy at $e^+e^-$ collider. (Right) Cross-section for the singly-charged quintuplet fermion production at $\gamma \gamma$ collider. }
\label{fig:eess2}
\end{center}
\end{figure}

An $e^+e^-$ linear collider can also be operated as a $\gamma \gamma$ and
an $e^- \gamma$ collider, as illustrated in Ref:\cite{osti_378863,DeRoeck:2003cjp}.
The highly intense photons for the collision are obtained by Compton back-scattering
laser photons on intense high-energy electron beams.
Due to the coupling with photon, charged particles can be produced with a considerably high cross-section at
these photon colliders. In the present model, the production of the
singly charged quintuplet fermions is possible via
$\gamma \gamma \rightarrow \Sigma^{+} \Sigma^{-}$,
$e^- \gamma \rightarrow \Sigma^{+} \phi{--}$ and
$e^- \gamma \rightarrow \Sigma^{-} \phi^{0^*}$ modes, along with the conjugate process in each case.
The production cross-section for the singly-charged quintuplet fermion at
$\gamma \gamma$ collider is shown in Fig:~\ref{fig:eess2} (right) as a function
of the center of mass energy.
These production modes, alone or combined with
$e^+e^-$ collision, show a great potential at future linear colliders. Over all,
the nonstandard decay modes of the quintuplet fermions offer different signals which
require alternate search strategies and there might be an opportunity for discovery and/or
exclusion at the HL-LHC and future linear colliders.

\vspace{-0.3cm}
\section*{Acknowledgments}
\vspace{-0.5cm}
N.K. acknowledges the support from the Dr. D. S. Kothari Postdoctoral scheme (201819-PH/18-19/0013). V.S. thanks Brajesh Choudhary for access to computers bought under the aegis of the Grant No. SR-MF/PS-0212014-DUB (G) of the DST (India).

\vspace{-0.3cm}
\bibliographystyle{unsrt}


text:
\section{Introduction}
  \input{intro.tex}

\section{Methods}
  \input{methods.tex}

\section{Evaluation}
  \input{results.tex}

\bibliographystyle{myrecomb}

\subsection{Background}
CoViT is based on ViT~\cite{dosovitskiy2020image}, a neural network primarily developed for image classification. We develop and apply a preprocessing step that employs MinHash~\cite{broder2000min}, whose role is extracting the informative feature vectors from the genome, which are further fed into a modified ViT for classification. In this section we provide a detailed introduction to the ViT and MinHash.

\subsubsection{Transformers}
\label{transformers}
Transformer is a deep neural network model proposed by Vaswani et al.~\cite{vaswani2017attention}. It became state of the art in many deep learning applications, such as  Natural Language Processing. More recently, the ViT which is an image classification deep neural network based on the transformer model, was introduced in ~\cite{dosovitskiy2020image}. ViT architecture is presented in Figure~\ref{fig:TIC}. It comprises an encoder and a multilayer perceptron (MLP) layer with a softmax activation function to perform predictions (in this paper we will refer to this MLP layer as MLP head). An encoder receives a sequence of representations, each of them is of dimensionality $d_{model}$, and outputs a sequence of learnt representations of the same dimensionality. For the sake of clarity we refer to the sequence of input representations as ``feature vectors sequence'' throughout this section. The transformer encoder layer is composed of two main sub-layers, namely the Multi Head Self Attention (MHSA) and the piece-wise MLP. In the following subsections we will present the structure and the purpose of those sub-layers.

\begin{figure*}[!t]
  \begin{center}
    \includegraphics[scale=0.3]{figs/TransformersInCoViT2.png}
  \end{center}
\caption{The Transformer encoder architecture is built as a stack of alternating $L$ multi-head self-attention (MHSA) layers and $L$ MLP blocks. Layer normalization is applied before every block and a residual connection is employed around each of the sub-layers (i.e., MHSA and MLP)}
\label{fig:TIC}
\end{figure*}

\textbf{Multi Head Self Attention:}
This sub-layer, depicted in Figure~\ref{fig:TIC} (right) receives as an input a sequence of feature vectors $f_1, f_2,\dots, f_n$, each of dimensionality $d_{model}$. The sub-layer outputs a sequence of learnt representations $r_1, r_2, \dots, r_n$, a representation per each feature vector. The output learnt representations are of the same dimensionality as the input ones.
The following presents the MHSA data flow. First, each feature vector $f_j$ is passed through $3h$ linear layers $\{q_l, k_l, v_l\}_{l=1}^{h}$, where $h$ denotes the number of heads, and 3 denotes the three types of different representations (i.e., query, key, and value). The purpose of those linear layers is to learn $h$ query representations for $f_j$, $h$ key representations for $f_j$ and $h$ value representations of $f_j$. Each linear layer defines a weight matrix as follows:
$W_l^q \in \mathbb{R}^{d_{model}\times d_{q}}$;\textbf{  } $ W_l^k \in \mathbb{R}^{d_{model}\times d_{k}}$;\textbf{  } $W_l^v \in \mathbb{R}^{d_{model}\times d_{v}}$.
They are used to calculate their respective representations:
$q_l(f_j) = f_j^T\cdot W_l^q$;\textbf{  } $k_l(f_j) = f_j^T\cdot W_l^k$;\textbf{  } $v_l(f_j) = f_j^T\cdot W_l^v$.

The query and key representations are used to calculate the degree of attention that the learnt output representation, $r_j$, will give to other feature vectors. The value representations are used to represent learnt properties of the feature vector that we intend to pass through the network. Formally, this first step of the MHSA layer is defined as follows:
\begin{align*}
    & A^l_j = \sum_{i=1}^{n}{similarity(q_l(f_j),k_l(f_i))}\cdot v_l(f_i)\\
    & similarity(q_l(f_j),k_l(f_i)) =
    softmax_i(\frac{q_l(f_j)^T k_l(f_1)}{\sqrt{d_k}}, \dots, \frac{q_l(f_j)^T k_l(f_1)}{\sqrt{d_k}})
\end{align*}
where $l$ refers to the $l$'th head, and $softmax_i$ is the $i$'th element of the softmax function, defined as
$softmax(x_1, x_2, \dots, x_n)=\frac{1}{\sum_{i=1}^{n}{e^{x_i}}}(e^{x_1}, e^{x_2},\dots, e^{x_n}).
$

Note that for a given head $l$, feature vectors with higher similarity to $f_j$ will have a larger dot product with $q_l(f_j)$ compared to other feature vectors. This drives the network to learn, for each representation, $r_j$, on which feature vectors of the input sequence it should focus its attention to learn the representations of value. To recap, representations $A_j^1, A_j^2 \dots, A_j^h$ depend on feature vectors that receive a high degree of attention.

Next, these representations are concatenated to a form $A_j = A_j^1 \circ A_j^2 \circ \dots \circ A_j^h$.

The last step of the MHSA is applying to $A_j$ a linear layer, to generate the MHSA output representation.

\textbf{Trainable weights and hyper-parameters:}
For each head, $l\in [1, h]$, we define three weight matrices $W_l^q, W_l^k, W_l^v$.
Additional weight matrix $W^o\in \mathbb{R}^{h\cdot d_v\times d_{model}}$ is reserved for the linear layer that processes $A_j$.

The MHSA hyper-parameters include the number of heads $h$, the dimensionality of the query and key representation $d_q=d_k$, and the dimensionality of the value representation $d_v$.

\textbf{Piece-wise Multi Layer Perceptron:}
The piece-wise MLP receives as an input a sequence of representations (feature vectors). For each feature vector it learns an output representation. Unlike the MHSA, in the piece-wise MLP layer each learnt representation depends only on its input feature vector and does not depend on other feature vectors in the sequence. This layer transforms each input feature vector into an intermediate representation of dimensionality $d_{ff}$, and learns the output representation which is of dimensionality $d_{model}$. Formally, to perform this calculation, we define two weight matrices
$
W^{ff}\in \mathbb{R}^{d_{model}\times d_{ff}}$ and \textbf{  } $W^m\in \mathbb{R}^{d_{ff}\times d_{model}}$.
and calculate the representation of $f$ as follows:
\begin{align*}
    MLP(f) = (ReLU(f^T\cdot W^{ff}))^T\cdot W^m
\end{align*}
\textbf{Trainable weights and hyper-parameters:}
There are two weight matrices $W^{ff}$ and $W^m$ as defined above. The only hyper-parameter is $d_{ff}$.

\textbf{Layer Normalization:}
Layer normalization~\cite{ba2016layer} is a method used to normalize the activities of the neurons of a layer improving the training speed for neural network models. It directly estimates the normalization statistics from the summed inputs to the neurons within a hidden layer.

\textbf{The Transformer Encoder:}
The Transformer Encoder (TE) comprises a stack of L identical encoder layers, each comprising an MHSA sub-layer followed by a piece-wise MLP sub-layer. Layer normalization is applied before each sub-layer, and a residual connection is employed around each of the two sub-layers.

\textbf{Hyper-parameters of the TE:}
 The number of encoder layers - $L$, the dimensionality of the query and key representations - $d_q = q_k$, the dimensionality of the value representation - $d_v$, and the dimensionality of the intermediate representation of the piece-wise MLP - $d_{ff}$.

\subsubsection{The MinHash Scheme}
The min-wise independent permutations (MinHash)~\cite{broder2000min} is a technique for similarity estimation.
MinHash was applied in many tasks in computational biology, including genome assembly~\cite{berlin2015assembling, shafin2020nanopore},
metagenomic gene clustering~\cite{rasheed2013map, muller2017metacache} and genomic distance estimation~\cite{ondov2016mash}. MinHash implements the following sketch function:

Given a set of characters $A$, a compression factor $n$, and a hash function $h$, the elements of the set $A$ are hashed using function $h$ to generate the set $H(A)$. Then the elements of $H(A)$ are sorted, and the smallest $n$ elements are returned as described in Algorithm~\ref{alg:sketch}.

\begin{algorithm}
\caption{The sketch algorithm $sketch(A)$}\label{alg:sketch}
\begin{algorithmic}
\Require set $A = \{a_1,\dots,a_{|A|}\}$, compression parameter $n$ and a hash function $h$
\State $sketch \gets \emptyset$
\State $H(A) \gets (h(a_1), \dots h(a_{|A|}))$
\State Sort H(A) to get $(h(a_{i_1}), h(a_{i_2}), \dots, h(a_{i_{|A|}}))$\\
\Return $(a_{i_1}, a_{i_2}, \dots, a_{i_{n}})$\footnotemark
\end{algorithmic}
\end{algorithm}
\footnotetext{These are the genome representative kmers referred to in Section 2.2.1}

The subset obtained by applying the sketch algorithm provides a good estimate for the Jaccard index~\cite{hancock2004jaccard}, defined as follows: Given two sets $A, B$ the Jaccard index is
$
J(A, B) = \frac{|A\cap B|}{|A\cup B|}
$.

Formally, given two sets, A and B, $\frac{|sketch(A) \cap sketch(B)|}{|sketch(A) \cup sketch(B)|} \approx J(A, B)$.
We use the sketch algorithm to extract feature vectors from a genome, as presented further.

\subsection{CoViT Architecture}

\begin{figure}[!t]
  \begin{center}
    \includegraphics[scale=0.48]{figs/archAndFEE.png}
  \end{center}
\caption{(a): A general overview of the CoViT Architecture.
(b): Feature Extractor workflow example.
}
\label{fig:covitarch}
\label{fig:fee}
\end{figure}

The top-level architecture of CoViT is shown in Figure~\ref{fig:covitarch}(a). The input for CoViT is an assembled genome in FASTA file format. To convert the genome into a sequence of feature vectors, we build a pipeline of preprocessing steps, comprising the feature extractor and the embedding layer. The feature extractor chooses a set of fragments (i.e., genome subsequences of fixed length), represents each fragment as a two-dimensional matrix and outputs those genome fragment matrices. Each genome fragment matrix is fed into an embedding layer consisting one neuron. This layer converts the genome fragment matrix into a genome fragment vector, designated the feature vector. The feature vectors are fed into the ViT which outputs the most likely placement candidates (by attaching probabilities to each lineage).

The feature extractor together with the embedding layer transform a genome into a sequence of representative feature vectors, which are the numerical representations of genome fragments. We implement the feature extractor using MinHash to preserve similarities in genome. It allows the feature extractor to find features in the query genome that are shared by the known lineages. The embedding layer is used as a dimensionality reduction step. It allows to reduce CoViT latency without affecting its placement accuracy.

\subsubsection{Feature extractor}
The first part of the CoViT pipeline receives the assembled genome and outputs matrices that represent genome fragments. Feature extractor employs MinHash to find similar fragments (i.e., features) in different genomes.

\begin{algorithm}
\caption{the Feature Extractor}\label{alg:featureExtractor}
\begin{algorithmic}
\Require genome $g$, compression parameter $n$, fragment size $f$, k-mer size $k$, hash function $h$
\State $G \gets \{ g_i \equiv g[i: i+k]\} \text{ } s.t.\text{ } i\in [0, n-k+1]$
\State $kmers \gets sketch(G, n, h)$
\State $left \gets floor(\frac{f-k}{2})$
\State $right \gets ceil(\frac{f-k}{2})$
\State $fragments \gets \{g[i-left, i+right] : g_i \in kmers\}$\\
\Return $one-hot-encode(fragments)$
\end{algorithmic}
{\footnotesize One hot encoding of a fragment refers to the operation of transforming each base in the fragment as follows:\newline $A \rightarrow [1, 0, 0, 0]$; \textbf{  }
$C \rightarrow [0, 1, 0, 0]$; \textbf{  }
$G \rightarrow [0, 0, 1, 0]$; \textbf{  }
$T \rightarrow [0, 0, 0, 1]$; \textbf{  }
$N \rightarrow [0, 0, 0, 0]$.\newline For the example please refer to Figure~\ref{fig:ele}(a).}
\end{algorithm}

An example is illustrated in Figure~\ref{fig:fee}(b). Input is a genome sequence $g$ of length $N=19$. We generate all possible kmers $G$ (where $k=4$), and extract ($n=3$) 4-mers from the genome using the MinHash scheme (i.e., $sketch(G,n=3,h)$). The next step is to generate fragments of length $f=8$ (by extending each of three 4-mers by $\frac{f-k}{2}=\frac{8-4}{2}=2$ bases in each direction). The last step of this workflow is to encode each basepair of the fragments using one-hot encoding.

For the genome sequence $g$ of length $N$, the extractor first generates a set $G$ of all possible kmers (genome sub-sequences of length $k$). It then sketches (i.e., applies MinHash \emph{sketch} function to) the set of kmers, $G$, to extract $n$ representative kmers of the genome (Algorithm \ref{alg:sketch}). Those kmers are used as anchors to be expanded to generate fragments. Last part is transforming a fragment into a numerical matrix where the genome basepairs A, G, C and T are encoded using one-hot encoding as described in Algorithm~\ref{alg:featureExtractor}.

\subsubsection{Embedding layer}
The input of the embedding layer is a genome fragment matrix of dimensions $f\times4$ (i.e., $f$ represents the fragment length and 4 is the dimensionality of the basepair one-hot encoding). It outputs a feature vector of size $f$. The embedding layer consists of one neuron which performs a base-wise linear transformation. Specifically, the embedding layer defines 5 learnable parameters, 4 weights $w_A, w_C, w_G, w_T$ and one bias term $b_N$. Given a genome fragment matrix, we transform each base (represented in one-hot encoding) to a numerical token as presented in Figure~\ref{fig:ele}(a), which  also contains an example in which we embed a fragment matrix of dimensions $8\times 4$.

\begin{figure}
  \begin{center}
    \includegraphics[scale=0.4]{figs/histAndELE.png}
  \end{center}
\caption{(a): Embedding a fragment matrix of dimensions $8\times4$. (b): The number of sequenced genomes per lineage histogram.
Most of the lineages (more than 800 out of 1,536 available) have less than 64 accessions available. The peak in the last bin accounts for all lineages with more than 960 accessions.
}
\label{fig:ele}
\label{fig:hist}
\end{figure}
\subsection{Setup}
In this section we describe the dataset, evaluation criteria, training process and the machine on which we train the model.

\subsubsection{Dataset}
We use the COVID-19 Data Portal~\cite{covid19data} as a source of sequenced SARS-CoV-2 genomes. There are 4,470,553 assembled genomes classified into 1,536 distinct lineages. The distribution of genomes across the lineages is highly non-uniform. For example, 463,426 samples are categorized as B.1.1.7, while the majority of the lineages has only a small number of samples, as depicted in Figure~\ref{fig:hist}(b). Such non-uniformity makes training the network challenging. To mitigate this, we select at most 1024 samples from each lineage and discard the rest when constructing the dataset. This way we prevent any potential over-fitting towards lineages like B.1.1.7. We also discard all samples that belong to lineages with fewer than 512 accessions, since such a small number of samples limits the model's ability to generalize. These steps form our dataset, which contains 240,734 samples belonging to 270 different lineage classes.

We divide our dataset into three non overlapping subsets: the training set, the validation set, and the test set. The training and validation sets are used to train and optimize the network. The test set is applied to evaluate CoViT's placement accuracy and to compare it with state-of-the-art solutions. Both the validation and test sets consist of 13,500 genomes randomly selected from the dataset. The rest of the samples comprise the training set.

\subsubsection{Classification accuracy and placement rate criteria}
We use the \emph{top-n accuracy} as CoViT placement accuracy criterion. It estimates the likelihood of the \emph{correct} placement (i.e., the lineage to which the query genome belongs to\footnote{Since the COVID-19 Data Portal assigns all sequenced genomes to their respective lineages, the correct results are known.}) appearing among the \emph{n} most probable results (i.e., \emph{n} lineages output by the network at the top of the descending-order probability list). Specifically, we measure the top-1, top-2 and top-5 accuracy of our model as a criteria of its performance. Such criteria make the comparison to the state-of-the-art straightforward.

By \emph{base ambiguity} we mean the ambiguous bases (i.e. R, Y, K, N etc.) in assembled SARS-CoV-2 genomes. Base ambiguity is a result of incomplete coverage. N marks the worst ambiguity because it can mean any base (A, G, C, or T). We further evaluate the accuracy of CoViT placement under different levels of base ambiguity. In our study, we conservatively use only the worst case (i.e., N) ambiguity.

Another criteria we apply in our analysis is the \emph{placement rate}, a fraction of accessions a placement algorithm is able to place (correctly or incorrectly) on the existing phylogenetic tree. The reasons for the placement rate to be lower than 100\% are base ambiguity and sequencing errors.

\subsubsection{Network Training}
We train CoViT using Adam~\cite{kingma2014adam} which is an optimized gradient descent technique, with $\beta_1 = 0.9, \beta_2=0.999$, batch size of 1024, and weight decay of 0.0001. Through various experiments, we found it most beneficial to choose the following hyper-parameters: $encoder repeats=4$, $d_{model}=256$, $d_{v} = 96$, $d_{k}= 96$, $d_{ff}=1536$, $h=18$. We also applied a dropout of rate 0.2 after each sub-layer.
The model was trained to minimize the categorical cross-entropy loss function. After fine-tuning and choosing those hyper-parameters, we train the model until the validation loss saturates.

\subsubsection{Hardware}
\label{HW}
We train CoViT on a desktop computer with Intel i7-9700K CPU with 8 cores, operating at 3.60GHz with 32GB of 2,666 MT/s DDR4 RAM, and NVIDIA's GeForce RTX 2080 Ti GPU with an 11GB / 14Gbps GDDR6 frame buffer, running at 1.545GHz.

\subsection{Results}

\subsubsection{Accuracy}
We measure CoViT performance using top-n accuracy.
We compare CoViT to the state-of-the-art phylogenetic placement tool, UShER~\cite{turakhia2021ultrafast}, which uses a maximum parsimony approach where it searches for a placement that requires the fewest additional mutations. UShER achieves 92.1\% top-1 accuracy over the test set, while CoViT presents a 2.1\% increment over UShER in top-1 accuracy performance, achieving 94.2\%. Moreover, CoViT achieves a top-2 accuracy of 97.9\% and top-5 accuracy of 99.8\% as summarized in  Table~\ref{Tab:t1}.

\begin{table}[hbt!]
\caption{Accuracy comparison
\label{Tab:t1}}
\begin{center}
{
\begin{tabular}{||c|c|c|c||}
 \hline
 Program & top-1 accuracy & top-2 accuracy & top-5 accuracy\\ [0.5ex]
 \hline\hline
 CoViT & 94.2\% & 97.9\% & 99.8\%\\ [1ex]
 \hline
 UShER & 92.1\% & - & -\\ [1ex]
 \hline
\end{tabular}
}
\end{center}
{\footnotesize
UShER uses another accuracy criterion called \emph{sister node placement}. They report the correct sister node placement in 97.2\% of cases. While this criterion resembles the top-2 accuracy, they are not identical, therefore we do not report it as a top-2 accuracy for UShER.}

\end{table}

\subsubsection{Run time and memory usage}
\label{efficiency}
We measure CoViT run time and memory usage and compare it to SARS-CoV-2 state-of-the-art placement tools. We compare CoViT and UShER, classifying the genomes in the test set, running on the machine specified in~\ref{HW} (see Table~\ref{Tab:t2}). We also present in Table~\ref{Tab:t2} the comparison to IQ-TREE multicore~\cite{nguyen2015iq}, and EPA-ng~\cite{barbera2019epa} running on a server with 160 processors (Intel Xeon CPU E7-8870 v.4, 2.10 gigahertz), each with 20 CPU cores. Those results where obtained while adding just one sequence to the SARS-CoV-2 global phylogeny containing 38,342 leaves~\cite{turakhia2021ultrafast}.

\newcommand{\mc}[2]{\multicolumn{#1}{c}{#2}}
\definecolor{Gray}{gray}{0.85}
\definecolor{LightCyan}{rgb}{0.88,1,1}

\begin{table}[hbt!]
\caption{Run-time and memory requirements comparison}
\label{Tab:t2}
\begin{center}
{
\begin{tabular}{|c|c|c|c|}
\hline
 Program & Avg run-time per sample & Peak memory used (GB) & speedup over UShER\\ [0.5ex]
 \hline
 \hline
 \rowcolor{LightCyan}CoViT (GPU) & \textbf{0.055s} & 5.02 & \textbf{3.53}\\
 \hline
 \rowcolor{LightCyan}CoViT (CPU) & 0.094s & 3.74 & 2.06\\
 \hline
 \rowcolor{LightCyan}UShER v0.5.6 & 0.195s & \textbf{0.49} & 1\\
 \hline
 \hline
 \rowcolor{yellow}IQ-TREE multicore v2.1.1 & 2791s & 12.85 & $7\cdot10^{-6}$\\
 \hline
  \rowcolor{yellow}EPA-ng v0.3.8 & 1658s & 791.82 & $6\cdot10^{-5}$\\
 \hline
\end{tabular}
}
\end{center}
{\footnotesize The results marked in light cyan are given running of the machine specified in section~\ref{HW}. The results marked in yellow are adopted from~\cite{turakhia2021ultrafast}.}
\end{table}

\subsubsection{Robustness to base ambiguity}
Some low quality SARS-CoV-2 genomes  may contain ambiguous bases, which is a result of incompletely covered genome or sequencing errors. It is highly desirable that such genomes are correctly placed even when the level of ambiguity is significant. To test the robustness of CoViT to base ambiguity, we mask  random sites (i.e., replace basepairs by N) in each SARS-CoV-2 genome of the test set, with the percentage of ambiguous bases ranging from 0 to 32\%. We then compare CoViT accuracy and placement rate on this \emph{ambiguous} test set (13,500 accessions with randomly masked sites) against the placement rate and accuracy of UShER. Figure~\ref{fig:ambig} presents the results.
At 32\% ambiguity, UShER's placement rate is 0.5\%, which means UShER places only 67 accessions out of 13,500. For the rest of accessions, UShER returns "failed to map" message.  For those accessions UShER is able to map, it reaches the top-1 accuracy of 60.5\%.
CoViT maintains a constant 100\% placement rate. At 32\% ambiguity, CoViT obtains the top-1 accuracy of 82.9\%.
CoViT maintains 100\% placement rate up to 50\% ambiguity (not shown in the figure), at which point CoViT's top-1 accuracy stands at 77.3\%. CoViT exhibits such resilience to base ambiguity because neural networks typically maintain high resilience when processing low quality (high error or ambiguity rates) data.

\begin{figure}[hbt!]
  \begin{center}
    \includegraphics[scale=0.32]{figs/ambig.png}
  \end{center}
\caption{Placement rate (left) and accuracy (right) as a function of the genome ambiguity. i.e., the fraction of the ambiguous bases in the input genome.
}
\label{fig:ambig}
\end{figure}

\section{Conclusion}
The massive amount of sequenced SARS-CoV-2 DNA material paves the way for real-time surveillance of the virus evolution and spread. At the same time, it puts a tremendous strain on phylogenetic analysis infrastructures. To prevent losing valuable genomic data and reaching incorrect conclusions, a new computing paradigm for phylogenetic analysis of viral pandemics is highly required.
CoViT is inspired by an enormous success of neural networks in implementing classification tasks.
CoViT extracts sets of features from a SARS-CoV-2 accession genome and uses Vision Transformer, an attention based neural network model, to accurately classify and place it in real-time. CoViT not only outperforms state-of-the-art tools both in accuracy and run-time, by placing a SARS-CoV-2 genome on a phylogenetic tree in 0.055s with 94.2\% accuracy. It also shows a strong resilience to low quality genome data (i.e. high ambiguity of the newly sequenced accessions). CoViT is able to place SARS-CoV-2 genomes with 50\% ambiguity (when N replaces A, G, C, and T in 50\% of the sites in each genome), while maintaining 77.3\% accuracy.

CoViT lifts the computational barrier to pandemic surveillance by improving the placement accuracy. It exhibits  highly robust performance under low quality and ambiguous genomic data while greatly decreasing the processing time. Thus it enables real-time evolutionary tracking of COVID-19 and future viral pandemics.

\subsubsection{Data set}
We use the COVID-19 Data Portal~\cite{covid19data} as a source of the sequenced SARS-CoV-2 genomes.
We sort the SARS-CoV-2 lineages based on the number of sequenced genomes (accessions) in each lineage. Neural network training requires a training set of a sufficient size to provide an accurate classification result. However, only a limited number of lineages have sufficient amount of sequenced genomes (accessions), thus limiting the number of lineages CoViT is able to accurately classify a newly sequenced SARS-CoV-2 genome to. The number of lineages with at least 256 accessions reaches only 375 out of 1536 recorded lineages. Most of the lineages has less than 50 accessions available as depicted in Figure~\ref{fig:hist}.

\subsubsection{Classification groups}
We create four groups of lineages based on the number of accessions in each, as follows:
\begin{itemize}
    \item CG1: 107 lineages, with at least 2048 sequenced genomes in each. To build the dataset for this group we take from each lineage exactly 2048 samples (i.e., sequenced genomes).
    \item CG2: 189 lineages, with at least 1024 sequenced genomes in each. To build the dataset for this group we take from each lineage exactly 1024 samples.
    \item CG3: 269 lineages, with at least 512 sequenced genomes in each. To build the dataset for this group we take from each lineage at least 512 samples, and 1024 samples at most.
    \item CG4: 375 lineages, with at least 256 sequence genomes in each. To build the dataset for this group we take from each lineage at least 256 samples, and 1024 samples at most.
\end{itemize}

The purpose of this division is to investigate how much the size of the training set affects the accuracy of the placement.

\subsubsection{Classification accuracy criterion}
CoViT classifies a newly sequenced (query) genome by sorting the known lineages in the descending order of their similarity to it. Therefore, one or two lineages of the highest similarity to the query genome are of the most interest.

We use the top-n accuracy as CoViT placement accuracy criterion. It estimates the likelihood of the \emph{correct} result (i.e. the lineage that fits the query genome) appearing among the \emph{n} highest-probability results (i.e. \emph{n} lineages placed by the network at the top of the descending-order similarity list). Specifically, we measure the top-1, top-2 and top-5 accuracy for every classification group.
Since the Covid-19 Data Portal assigns all sequenced genomes to their respective lineages, we know during the network training which results are correct, and therefore are able to tag them.

\subsubsection{Network training}
\label{Network training}
We train CoViT for each classification group, using the following training and validation datasets:
\begin{itemize}
    \item 80\% of the sequenced genomes are applied as the training set
    \item Remaining 20\% are employed as the validation set
\end{itemize}
We train CoViT on NVIDIA's GeForce RTX 2080 Ti GPU with an 11GB / 14Gbps GDDR6 frame buffer, running at 1.545GHz. CoViT was implemented using the TensorFlow framework in Python. To train the network, we use two main techniques, the Greedy Layer-Wise Pretraining~\cite{bengio2006greedy} and Transfer Learning~\cite{torrey2010transfer}.

\paragraph{Greedy Layer-Wise Pretraining}
Training deep neural networks is challenging due to the vanishing gradient problem (i.e. the weights in layers close to the input layer are not updated in response to errors calculated on the training dataset). The Greedy Layer-Wise Pretraining allows deep neural networks to be successfully trained, achieving state-of-the-art performance. Pretraining involves adding a new hidden layer to a model and refitting, allowing the new model to learn the inputs from the existing hidden layers, while keeping the weights for the existing hidden layers fixed and not trainable. This technique is based on the observation that training a shallow network is easier than training a deep one~\cite{bengio2006greedy}.

\paragraph{Transfer Learning}
A technique that focuses on storing knowledge gained while solving one problem, and applying it to a different but related problem~\cite{west2007spring}. In transfer learning, we first train a base network on a base dataset and task, and then we re-purpose the learned features, or transfer them, to a second target network to be trained on a target dataset and task. This works well if the features are general, meaning suitable to both base and target tasks, rather than being specific to the base task~\cite{yosinski2014transferable}.

\paragraph{Regularization}
To generalize the results on the training set better, we use several regularization techniques in parallel. We employ a dropout after each MSA and MLP sub-layers in the network. Moreover we use L2 regularization~\cite{krogh1991simple} on the trainable weights of the network (not including the biases).

\subsection{Evaluation results}
\subsubsection{Classification accuracy}
We started by building a CoViT model for the CG3 classification group. To develop deeper networks and achieve best results we chose the following hyper-parameters
\begin{center}
\begin{tabular}{||c|c|c|c|c|c|c||}
 \hline
 $d_{model}$ & $d_{v}$ & $d_{k}$ & $d_{ff}$ & $heads$ & $dropout\_rate$ & $regularization\_coef$ \\ [0.5ex]
 \hline\hline
 256 & 96 & 96 & 1536 & 18 & 0.1 & 0.0001 \\ [1ex]
 \hline
\end{tabular}
\end{center}
Also, we chose the $kmer$ length to be $16$, the fragment length to be $256$ (i.e., equal to the $d_{model}$) and the comperession factor, $n$, is set at $256$. These parameters are applied to all trained models.
We started training a model for the CG3 classification group using Supervised Greedy Layer-Wise Pretraining~\cite{bengio2006greedy}. We started the pretraining process by building a base model consisting of two TE layers. We trained the model for 60 epochs with a batch size of 256 \emph{examples} (i.e., assembled genomes). The average training time per epoch was 8.58[min], achieving the top-1 accuracy of 79.3\% on the validation set. Then we added two new TE layers to the base model, generating a 4-layer TE model. This 4-layer model was trained for 120 epochs, with the batch size of 256 examples. The average training time per epoch was 16.84[min]. We achieved the top-1 accuracy of 92.7\%, top-2 accuracy of 97.2\% and top-5 accuracy of 99.1\% on the validation set. The accuracy as a function of the number of epochs is depicted in Figure~\ref{fig:accs}(c).

To build the CoViT model for the CG1, CG2 and CG4 classification groups we used transfer learning~\cite{torrey2010transfer}. More specifically, we took the 4-layer learnt model for the CG3 as the base for the CG2 and the CG4 models, where the only component to change was the MLP head (presented in section~\ref{transformers}). We then trained each model for a number of epochs specified in Table~\ref{Tab:t1}, with a batch size of 256 examples.

The accuracy as a function of the number of epochs for these remaining classification groups is presented in Figure~\ref{fig:accs}(a), (b) and (d).
Following are the validation set results. The model for CG1 classification group achieves the highest top-1 accuracy of 97.1\%, top-2 accuracy of 99.0\% and top-5 accuracy of 99.8\%. Respectively, the model for CG2 classification group achieves the highest top-1 accuracy of 95.8\%, top-2 accuracy of 98.1\% and top-5 accuracy of 99.8\%.
The model for CG3 classification group achieves the highest top-1 accuracy of 92.7\%, top-2 accuracy of 97.2\% and top-5 accuracy of 99.1\%. The model for CG4 classification group achieves the highest top-1 accuracy of 92.1\%, top-2 accuracy of 96.9\% and top-5 accuracy of 99.1\%.
These accuracy results and additional training information are summarized in Table~\ref{Tab:t1}.

\begin{figure}[!t]
  \begin{center}
    \includegraphics[scale=0.43]{figs/Accuracy_fig6.png}
  \end{center}
\caption{Accuracy vs. number of epochs for different models and four classification groups: CG1 (a), CG2 (b), CG3 (c) and CG4 (d).
}
\label{fig:accs}
\end{figure}

\begin{table}[ht]
\caption{Accuracy results and training information of the different models on the different classification groups\label{Tab:t1}}
\begin{center}
\begin{tabular}{||c|c|c|c|c|c|c||}
\hline
classification group & transferred from & epochs & avg $\frac{train\_time}{epoch}$ & top-1 acc & top-2 acc & top-5 acc\\ [0.5ex]
\hline\hline
CG1 & CG2 & 15 & 12.89min & 97.1\% & 99.0\% & 99.8\% \\
\hline
CG2 & CG3 & 40 & 13.09min & 95.8\% & 98.1\% & 99.8\% \\
\hline
CG3 & - & 60, 120 & 8.58min, 16.84min & 92.7\% & 97.2\% & 99.1\% \\
\hline
CG4 & CG3 & 155 & 21.52min & 92.1\% & 96.9\% & 99.1\% \\ [1ex]
\hline
\end{tabular}
\end{center}
\end{table}

\subsubsection{Training and inference times}
The following Table~\ref{Tab:t2} summarizes the training and inference times for each classification group. Inference time for the CG1 is 1.45[ms] which is almost three orders of magnitude faster than probably the fastest state of the art placement tool UShER~\cite{turakhia2021ultrafast}.

\begin{table}[ht]
\caption{Training and inference time of the various models.\label{Tab:t2}}
\begin{center}
\begin{tabular}{||c|c|c||}
 \hline
 classification group & train time per example & inference time per example\\ [0.5ex]
 \hline\hline
 CG1 & 4.59ms & 1.45ms\\
 \hline
 CG2 & 4.67ms & 1.47ms\\
 \hline
 CG3 & 4.70ms & 1.49ms\\
 \hline
 CG4 & 4.71ms & 1.49ms\\ [1ex]
 \hline
\end{tabular}
\end{center}
\end{table}

\subsubsection{Comparison with state-of-the-art}
Several methods have been developed for the purpose of placing genome sequence samples onto an existing reference phylogeny tree~\cite{minh2020iq, barbera2019epa, loytynoja2012accurate, ruan2007treefam}. However, these methods are far too slow to enable real-time genomic placement.
To address the speed of interpretation of viral transmission and evolutionary dynamics, ultrafast bootstrapping~\cite{hoang2018ufboot2, minh2013ultrafast} have been recently developed. While those methods provide significant speedup, they are not applicable to the problem of placing individual samples onto a reference phylogeny.

Below we address several state of the art classification and placement tools, including Kraken~\cite{wood2019improved}, PACIFIC~\cite{acera2021pacific}, UShER~\cite{turakhia2021ultrafast} and PangoLearn~\cite{o2021assignment}.

Kraken is one of the most popular genome classifiers. It uses exact matching of kmers of the genome(s) in question against a known organisms' genomes database. Kraken was not specifically designed to classify SARS-CoV-2 (or another viral pathogen with pandemic potential) lineages and therefore is not sensitive enough to place SARS-CoV-2 variants which differ from one another by only several basepairs.

PACIFIC is arguably the first deep classification neural network based genome classifier. PACIFIC has been shown to classify genomes into a very limited group of classes (5 distinct virus classes and a human genome class).

PangoLearn is an algorithm for lineage assignment which relies on machine learning and uses random forests. The model was trained using SARS-CoV-2 sequences from GISAID~\cite{gisaid}.

\begin{figure}[!t]
  \begin{center}
    \includegraphics[scale=0.25]{figs/covitvsusher.png}
  \end{center}
\caption{CoViT accuracy of the different models against UShER on a randomly selected set of accessions.
}
\label{fig:cvu}
\end{figure}

We compare the top-1 accuracy and placement speed of four CoViT models  (i.e., CG1, CG2, CG3, CG4) against UShER~\cite{turakhia2021ultrafast} by placing 13,500 genome assemblies randomly selected from the Covid-19 Data Portal~\cite{covid19data}.
UShER is one of the latest phylogenetic placement tools. It uses a maximum parsimony approach where it searches the entire reference variant-annotated tree
for a placement that requires the fewest additional mutations.
One reported limitation of UShER is that it only supports mutations caused by replacements (does not support mutations caused by indels).

For models CG1 through CG4, CoViT achieves the average placement accuracy of 96.1\%, 94.2\%, 93.8\% 93.9\% respectively, while UshER's average accuracy is 92.1\% as depicted in Figure~\ref{fig:cvu}. To test the placement speed, we ran both CoViT and UShER on Intel i7-9700K CPU with 8 cores, operating at 3.60GHz. Additionally, CoViT was implemented on NVIDIA's GeForce RTX
2080 Ti GPU. The placement times for these 13,500 genome assemblies are presented in table~\ref{Tab:t3}. The rest of the table is adopted from \cite{turakhia2021ultrafast}.

\begin{table}[t]
\caption{Placement time of 13,500 SARS-CoV-2 genomes assemblies.
\label{Tab:t3}}
\begin{center}
\begin{tabular}{||c|c|c||}
 \hline
 program & time & speedup vs UshER\\ [0.5ex]
 \hline\hline
 CG1 (GPU) & 12.42m & 3.6\\
 \hline
 CG2 (GPU) & 12.45m & 3.6\\
 \hline
 CG3 (GPU) & 12.47m & 3.6\\
 \hline
 CG4 (GPU) & 12.52m & 3.5\\
 \hline
 CG1 (CPU) & 21.35m & 2.1\\
 \hline
 CG2 (CPU) & 21.35m & 2.1\\
 \hline
 CG3 (CPU) & 21.35m & 2.1\\
 \hline
 CG4 (CPU) & 21.38m & 2.1\\
 \hline
 UshER & 44.1m & 1\\ [1ex]
 \hline
\end{tabular}
\end{center}
\end{table}

redpajama_set_name: RedPajamaArXiv