# Study of the $tH$ production at 14TeV and 100TeV $pp$ colliders

Yi-Dong Song<EMAIL_ADDRESS>School of Physics, Shandong University,
Jinan, 250100, P. R. China  Shou-Shan Bao<EMAIL_ADDRESS>Institute of
Frontier and Interdisciplinary Science, Shandong University, Qingdao, 266237,
P. R. China

###### Abstract

The flavor-changing neutral current interactions in the standard model are
suppressed seriously and such interactions can be used to search the new
physics beyond SM. The top quark and Higgs bosons are heavier than the other
particles in SM, we can expect the new physics plays a more important role in
their interactions. In this work, we study the flavor- changing neutral
current interactions between the top quark $\bar{t}qH$ through the production
of the single top associated with a Higgs boson on 14TeV and 100TeV $pp$
colliders. We consider the leptonic decay channels of Higgs and study the
signal. We find a sensitive region $0.4\leq\Delta R\leq 1.4$ for the leptons
from Higgs bosons and $\Delta R\geq 1.8$ between the jets from top quark and
leptons. We investigate the detective abilities of the hadron colliders for
the processes. Although this process seems not comparable with $pp\to
t\bar{t}$ with $t\to Hq$ decay, it is still attractive since this process can
be used to distinguish the $y_{tu}$ and $y_{tc}$.

##  I Introduction

The discovery of the Higgs boson at the Large Hadron Collider (LHC)[1, 2] is a
great triumph of the Standard Model (SM). Following that, one can expect to
test and search for new physics through the precision measurement of the Higgs
boson properties. In SM, the Higgs doublet is introduced to generate masses
for the gauge bosons and the fermions. Most parameters in SM are related to
Yukawa couplings. Since the top quark and Higgs boson are heavier than the
other particles in SM, we can expect that high energy new physics may play
more important roles in their couplings and the Yukawa coupling between the
top quark and Higgs boson could provide better clues to study the mechanism of
electroweak symmetry breaking and new physics beyond SM. On the other hand, so
far there is no signal of new physics particles has been found in direct
searches. Therefore indirect search through the precision measurement can be
expected to constrain the new physics scale or coupling. In SM the flavor-
changing neutral current (FCNC) couplings are absent at the tree level and
suppressed at the loop level. So that such processes are thought to be good
clues to search for the new physics beyond SM indirectly.

The FCNC interactions of top quark are small and the branch ratios of such
decay are predicted as $Br(t\to Hc)\sim 10^{-14}$ and $Br(t\to Hu)\sim
10^{-17}$ in SM[3, 4], which are far below the sensitivity of LHC. Many new
physics extensions to the SM can enhance such interaction at tree level or
loop level, e.g. two Higgs doublet model[5, 6, 7], supersymmetry model[8, 9,
10] and other new physics models[11]. One can expect that the new physics
contributions would play significant roles in such FCNC interactions and such
contributions could be discovered or constrained on high energy colliders.

In this work, we focus on the effective FCNH interaction between the Higgs
boson and the top quark which can be expressed in general as [14, 15]

$-\mathcal{L}_{Y}\supset\sum_{q=u,c}y_{tq}\bar{q}tH+h.c.,$ (1)

where the $y_{tq}$ is the FCNC Yukawa coupling constant. Such anomalous
couplings have been studied widely. Many works have been done through $t\to
Hq$ rare decays [15, 16, 17], the single top quark production with Higgs boson
[18, 19, 20], and same sign top quark production [21]. Recently, a study
invoving $pp\to tH$ and $pp\to tqH$ on 100TeV hadron collider is done in Ref.
[22]. Through $pp\to t\bar{t}$ with $t\to Hq$ channel, ATLAS[23] and CMS[24]
set upper limits at the 95% confidence level (C.L.) on the couplings as

ATLAS: $\displaystyle Br(t\to Hc)<9.9\times 10^{-4},$ (2) $\displaystyle
Br(t\to Hu)<7.2\times 10^{-4},$ (3) $\displaystyle\rm{CMS:\leavevmode\nobreak\
}$ $\displaystyle Br(t\to Hc)<9.4\times 10^{-4},$ (4) $\displaystyle Br(t\to
Hu)<7.9\times 10^{-4}.$ (5)

Since the $u$-quark and $c$-quark from top quark decay can be distinguished,
the LHC can not distinguish the vertexes $\bar{t}uH$ and $\bar{t}cH$. However,
the two vertexes can be distinguished in $tH$ and $\bar{t}H$ productions,
since the PDF luminosity of the $u$-quark in the proton is larger than the sea
quark.

In this work we focus on the process $pp\to t/\bar{t}H$ at high luminosity LHC
and 100TeV $pp$ colliders in future [13, 13]. In Ref.[19], the QCD NLO
corrections have been studied at LHC. In Ref. [12] this process was analyzed
with $H\to b\bar{b}$ channel, and it was found that the angular of the final
states could be used to separate the boosted Higgs signal from the main
background. Since the isolated high $p_{t}$ electrons and muons are typically
clean and easy detected on colliders, we focus on the $h\to VV^{*}$ ($V=W,Z$)
channels to study the $tH$ associated production at hadron colliders and try
to examine whether such channels are prossiming channels for measurement of
the anomalous Yukawa couplings at $pp$ colliders.

This work is organized as follows. In Sec. II, the $H\to
WW^{*}\to\ell^{+}\ell^{-}+E\\!\\!\\!/_{T}$ channel is analyzed for 14TeV and
100TeV $pp$ colliders. In Sec. III, the $H\to ZZ^{*}\to 4\ell^{\pm}$ channel
is argued. The results and the discussion are given in Sec. IV and a short
summary is given in Sec. V.

##  II $h\to WW^{*}\to 2\ell^{\pm}+E\\!\\!\\!/_{T}$

We study the anomalous interaction in Eq. (1) through $pp\to t/\bar{t}H$ and
the Feynman diagrams are shown in Fig. 1. Such a process can also be induced
by the FCNC coupling of strong interaction $\bar{t}qg$ in some new physics
models. However, there will be only $s$-channel contribution which is
suppressed at high energy colliders.

Figure 1: The Feynman diagram for the $tH$ and $\bar{t}H$ production at hadron
colliders.

In this work, we focus on the leptonic decay channels of Higgs and hadronic
decays of the top quark. Firstly we consider the $h\to WW^{*}$ channel. The
final states are characterized by two charged leptons and three jets including
one $b$-jet and two light jets. The dominate background is the process $pp\to
t\bar{t}j\to WWb\bar{b}j$ which is denoted as $t\bar{t}$ in this work. Also
the processes of EW gauge bosons produced with light jets $W^{+}W^{-}+3j$ and
$pp\to Z+3j$ are considered, which can be suppressed by the b-tag and the
invariant mass of the two charged leptons.

We generate the events of the signals and backgrounds with Madgraph [25] with
Pythia8 [26] and Deltphes3 [27]. The default values in Madgraph for the SM
parameters are used. The build-in ATLAS card in Delphes is used for the
simulation of 14 TeV hadron collider and the FCC-h card is used for 100TeV
hadron collider. We choose NNPDF23LO [28] parton distribution function set in
this work.

In the analysis we need $N_{j}\geq 3,N_{l}\geq 2$ and $N_{b}=1$, and some
other basic acceptance cuts are needed for the selection of events,

$\displaystyle P_{T}\geq 25\mathrm{GeV},\quad E\\!\\!\\!/_{T}>25\mathrm{GeV},$
(6a) $\displaystyle\Delta R\geq 0.4,\quad\lvert\eta\rvert\leq 2.5,$ (6b)

where the $p_{T}$ and $\eta$ are the transverse momentum and rapidity of
leptons and jets. $E\\!\\!\\!/_{T}$ is the missing transverse energy. The
angular distance $\Delta R=\sqrt{(\Delta\eta)^{2}+(\Delta\phi)^{2}}$ is the
separation in $\eta-\phi$ plane between different final states. Indeed, the
acceptance cuts depend on the detector could be different for the jets and
leptons. Here, we employ the same cuts for different final states for
simplicity. All the above cuts are denoted as Cut-I in our analysis.

To reconstruct the top quark and Higgs boson, we need other selection
conditions. In Fig. 2, the distributions for the invariant mass of the three
jets and that of two charged leptons are shown. Since the $b$-jet and two
light jets are from one top quark, we can expect the invariant mass $m_{bjj}$
to be close to the top mass. Furthermore, the two light jets are from one $W$
boson, so we can employ a cut on the invariant mass $m_{jj}$ near $M_{W}$. The
three jets in the $t\bar{t}j$ background are from two different top quarks
where one of the two $b$-jets is miss tagged. Since then there is no peak in
the distributions of the backgrounds. According to the plots, we can apply the
cuts for the $m_{bjj}$ and the $m_{jj}$ to reconstruct the top quark as

$\displaystyle 130\mathrm{GeV}\leq m_{bjj}\leq 200\mathrm{GeV}.$ (7a)
$\displaystyle 50\mathrm{GeV}\leq m_{jj}\leq 100\mathrm{GeV}.$ (7b)

The two charged leptons which are from the Higgs boson, carry a part of the
Higgs energy while the other part is carried by the neutrinos. Since then
there will be $m_{\ell\ell}<M_{W}$. However, there will be a narrow peak near
$M_{Z}$ in the $m_{\ell\ell}$ distribution for the background $Zjjj$. It’s
easy to veto such events with an up limit cut on the $m_{\ell\ell}$ as

$\displaystyle m_{\ell\ell}\leq 81\mathrm{GeV}.$ (8)

We categorize the cuts of $m_{bjj}$, $m_{jj}$ and $m_{\ell\ell}$ as Cut-II.

Figure 2: The normalized distributions of the invariant mass (GeV) of the
three jets, two light jets and two charged leptons at 14TeV and 100 TeV $pp$
colliders.

The opening angle between the two charged leptons tends to be small due to the
$V-A$ structure of the interaction between the leptons and $W$ boson [30].
This can be used to separate the signal from the backgrounds where the charged
leptons are likely to have a large opening angle. Furthermore, the Higgs and
the top quark are boosted strongly if the collision energy is high. In the
situation, the three jets are also close to each other and the leptons from
Higgs should be moving more closely. In the center-of-mass frame, the top
quark and Higgs boson are moving back to back, therefore the angular distance
between the leptons and the jets would be far. In Fig. 3, we show the
distributions for $\Delta R_{\ell\ell}$ and $\Delta R_{\ell j}$. According to
that, we employ the following cuts which are denoted as Cut-III,

$\displaystyle\Delta R_{\ell\ell}\leq 1.4,\quad\Delta R_{\ell j}\geq 1.8.$ (9)

Figure 3: The normalized distributions of $\Delta R_{\ell\ell}$ and $\Delta R_{\ell j}$ at 14 TeV and 100 TeV $pp$ colliders.  | Parton level | Cut-I | Cut-II | Cut-III | $\epsilon_{III}$
---|---|---|---|---|---
$gc$ | 72.37fb | 4.652fb | 0.7024fb | 200.7ab | 28.6%
$gu$ | 504.4fb | 25.62fb | 3.306fb | 1241ab | 37.6%
$t\bar{t}j$ | 56.23pb | 2866fb | 43.97fb | 6148ab | 14.0%
wwjjj | 2.427pb | 14.24fb | 0.08803fb | 9.266ab | 10.6%
zjjj | 687.5pb | 992.6fb | 9.499fb | - | -
Table 1: The cross-sections of the signals and backgrounds with cuts on 14TeV $pp$ collider. The $gc$ ($gu$) denote the signal processes $gc\to tH$ and $g\bar{c}\to\bar{t}H$ ($gu\to tH$ and $g\bar{u}\to\bar{t}H$). The acceptance of the Cut-III $\epsilon_{III}$ is listed in the last column. The ”-” means there is no event left after employing the cut, which requires more events to be generated for the analysis. | Parton level | Cut-I | Cut-II | Cut-III | $\epsilon_{III}$
---|---|---|---|---|---
gc | 2.711pb | 0.1554pb | 22.16fb | 8.839fb | 39.9%
gu | 9.115pb | 0.2354pb | 36.34fb | 16.65fb | 45.0%
$t\bar{t}j$ | 626.8pb | 14.89pb | 127.9fb | 13.79fb | 10.8%
Table 2: The cross-sections of the signals and the main background with cuts
on 100TeV hadron collider.The convention of the symbols is the same as that in
Table. 1.

The cross-sections for signals and backgrounds after the cuts are given in
Table. 1 for 14 TeV and Table. 2 for 100 TeV $pp$ colliders. Only the
dominated background $t\bar{t}j$ is represented for 100TeV collisions. All the
events are generated at the tree level. And to avoid the soft-collinear
divergence, we employ some loose parton-level cuts on the final quarks and
leptons in Madgraph. These cuts are

$\displaystyle p^{i}_{T}\geq 10\mathrm{GeV},\quad\Delta R_{ij}\geq 0.2.$ (10)

where $i,j$ denote all the leptons and quarks in the final states. The ”-”
means there is event left after employing the cuts. To obtain the cross-
section, more events should be generated. Indeed we generate 290k events for
the background $Zjjj$. Since there is no event left with the cuts, we can
estimate a upper limit for the cross-section as 2.375fb with the event number
$\leq 1$. The acceptance for Cut-III is defined as the ratio between the
cross-sections with and without the Cut-III. From that we can see that the
cuts of the $\Delta R$s work better on 100 TeV than that on LHC. That can be
understood as more Higgs bosons tends to be boosted at higher energy $pp$
collisions.

##  III $h\to ZZ^{*}\to 4\ell^{\pm}$

In this section, we make an analysis of the four charged leptons final states
of the Higgs boson. The branching ratio of the Higgs decay $H\to ZZ^{*}\to
4\ell$ is 2 orders smaller than that of $H\to WW^{*}\to 2\ell$. This makes
that such a channel is not promising on LHC. The final states are
characterized by four charged leptons, one b-jet and two light jets. Such
final states can be from $ZZ+X$ production in SM. However, the b-tag will
suppress this background strongly. Here the background we consider are $pp\to
ZZ+b\bar{b}+X$ and $pp\to ZV^{*}+b\bar{b}+X$ where $V=\gamma,Z$. We denote the
total background as $Zllbbj$. Similar to the analysis in the above section, we
need $N_{\ell}\geq 4$, $N_{j}\geq 3$ and $N_{b}=1$ and the same requirement as
Eq. (6). These cuts are denoted as Cut-I.

Moreover, we need the cut on the jets invariant mass $m_{3j}$ and the leptons
invariant mass $m_{4l}$ to reconstruct the top quark and the Higgs boson. The
distributions for the signals and backgrounds are shown in Fig. 4. Similar to
the previous analysis, we choose a cut for the invariant mass of the two light
$M_{jj}$ near $M_{W}$ and a cut for the invariant mass of the three jets
$m_{bjj}$ near the top quark mass. There are two pairs of charged leptons in
the signal, and one of them are from an on-shell $Z-$boson. We select the pair
of charged leptons with same flavor whose invariant mass is most close to
$M_{Z}$ to be the leptons from the on-shell $Z-$boson.

$\displaystyle
m^{r}_{z}=\\{m_{\ell^{+}\ell^{-}}|\mathrm{min}(|m_{\ell^{+}\ell^{-}}-M_{Z}|)\\}.$
(11)

We collect the cuts as follows and named as Cut-II,

$\displaystyle 130\mathrm{GeV}\leq m_{3j}\leq 200\mathrm{GeV}.$ (12a)
$\displaystyle 50\mathrm{GeV}\leq m_{jj}\leq 100\mathrm{GeV}.$ (12b)
$\displaystyle 110\mathrm{GeV}\leq m_{4\ell}\leq 140\mathrm{GeV}.$ (12c)
$\displaystyle 80\mathrm{GeV}\leq m^{r}_{z}\leq 100\mathrm{GeV}.$ (12d)

Figure 4: The normalized distributions for signal and background at 14TeV and 100TeV hadron colliders. | Parton level | Cut-I | Cut-II | Cut-III
---|---|---|---|---
14TeV
gc | 1.554fb | 2.694ab | 0.5284ab | 0.1036ab
gu | 10.80fb | 13.03ab | 2.303ab | 0.5758ab
$Z\ell\ell b\bar{b}j$ | 21.30ab | 5.705zb | 129.7yb |  $\sim$0.3007yb
100TeV
gc | 58.12fb | 111.2ab | 23.63ab | 3.100ab
gu | 195.4fb | 216.3ab | 49.51ab | 10.42ab
z$\ell\ell b\bar{b}$j | 48.1fb | 509.0ab |  $\sim$208.3zb |  $\sim$42.66yb
Table 3: The cross-sections of the signals and backgrounds after employing the
cuts. For efficiency, we employ $100\mathrm{GeV}\leq m_{4\ell}\leq
150\mathrm{GeV}$ cut at parton level in Madgraph in addition to Eq. (13). The
backgrounds are estimated roughly as explained in the text.

The four lepton in the background $ZZ+jets$ are from two on-shell $Z$ boson,
therefore their invariant mass is greater than $2M_{Z}$ as shown in Fig. 4.
The Cut-II can veto most events of the process $ZZ+jets$. Finally, we employ
Cut-III which is the same as Eq. (9). The cross-sections of the signals and
backgrounds are shown in Table. 3. One can see that the cross-sections of the
signals with cuts seem too small to be detected on LHC. However, this channel
could be detectable on 100TeV $pp$ colliders with high luminosity. It should
be clarified that the cross-sections for the backgrounds here are estimated
very roughly. We generated 328476 events for LHC and 38258 events for 100TeV
colliders with some parton-level cuts which are given as follows

$\displaystyle p^{i}_{T}\geq 10\mathrm{GeV},\quad|\eta^{i}|\leq 5,\quad\Delta
R_{ij}\geq 0.25.$ (13)

And there is no event left with Cut I+II. To estimate the background with all
the cuts, firstly we obtain the acceptance rates of three cuts separately and
then assume the final acceptance rate as the product of the three cut
acceptance rates.

##  IV Results and discussion

From the previous analysis, the $WW^{*}$ channel is a promising channel to
detect the FCNH couplings. Now we try to illustrate the detective potential of
LHC and the hadron collider in the future. The likelihood function is defined
as the Poisson probability [31]

$\displaystyle L(n,n_{b}+n_{s})=\frac{(n_{b}+n_{s})^{n}}{n!}e^{-n_{b}-n_{s}}$
(14)

where the $n$, $n_{b}$ and $n_{s}$ denote the observed data, SM background and
signal events. In doing so we neglect the systematic uncertainties of the
signal and the background estimation. To get the expected upper limit, we can
replace $n$ with $n_{b}$. The significance of the signal events relative to
the background events is defined as

$\displaystyle\sigma_{s}(y_{tq})=\sqrt{-2\ln\frac{L(n_{b},n_{b}+n_{s})}{L(n_{b},n_{b})}}\simeq\frac{n_{s}}{\sqrt{n_{s}+n_{b}}}.$
(15)

The values of $y_{tq}$ corresponding to $\sigma_{s}\leq 2(3)$ can be regarded
as the excluded region at 95%(99%) confidence level. As shown in Table. 1,
Table. 2 and Table. 3, the $WW^{*}$ channel is more guaranteed than the
$ZZ^{*}$ channel, and we can combine the $WW^{*}$ and $ZZ^{*}$ channels with
the product of the likelihoods defined in Eq. (14).The plots of the upper
limits of the $y_{tu}$ and $y_{tc}$ depend on the luminosity of the hadron
collider are given in Fig. 5 and Fig. 6.

Figure 5: The upper limits of $y_{tu}$ and $y_{tc}$ depend on the luminosity
at LHC.

Figure 6: The upper limits of $y_{tu}$ and $y_{tc}$ depend on the luminosity
at 100TeV $pp$ collider.

Combing the $WW^{*}$ and $ZZ^{*}$ channels, the upper limits of the FCNH
couplings on $pp$ collider with high luminosity $\mathcal{L}=3000/fb$ at 95%
C.L. can be obtained as

14TeV: $\displaystyle y^{2}_{tu}\leq 6.96\times 10^{-2},\quad y^{2}_{tc}\leq
0.429,$ (16) 100TeV: $\displaystyle y^{2}_{tu}\leq 8.18\times 10^{-3},\quad
y^{2}_{tc}\leq 0.0154.$ (17)

We can see that this channel at LHC is not as good as the $H\to b\bar{b}$
channel studied in Ref. [12]. Since then we combine all the three channels of
$H\to b\bar{b}$, $WW^{*}$ and $ZZ^{*}$ and obtain the upper limits on LHC for
$\mathcal{L}=3000/fb$ as,

$\displaystyle y^{2}_{tu}\leq 1.08\times 10^{-3},\quad y^{2}_{tc}\leq
7.08\times 10^{-3},$ (18)

These constraints can be represented as the upper limits of the top decay,

$\displaystyle\mathrm{Br}(t\to Hc)\leq$ $\displaystyle 4.09\times 10^{-3},$
(19) $\displaystyle\mathrm{Br}(t\to Hu)\leq$ $\displaystyle 6.06\times
10^{-4}.$ (20)

According to the study in Ref. [32], the constraints on the FCNC decays of the
top quark could be reached to $Br(t\to Hu)\leq 1.2\times 10^{-4}$ and $Br(t\to
Hc)\leq 1.0\times 10^{-4}$ at 95% C.L. on the LHC with $\mathcal{L}=3000/fb$.
From that, we can see that the process $pp\to t/\bar{t}H$ seems not
completable with the direct search for the decay $t\to Hq$ on LHC. However,
the couplings $y_{tu}$ and $y_{tc}$ can not be well distinguished in the top
decay on the hadron colliders. Once the FCNC decay mode $t\to Hq$ is measured,
the $pp\to t/\bar{t}H$ can be used to identify the interaction of $y_{tq}$.

##  V Summary

In this work we studied the FCNC Yukawa couplings of the top quark at 14TeV
and 100TeV hadron colliders through the process $pp\to t/\bar{t}H$. We analyze
the leptonic decay channels of the Higgs boson. We find that the special
region where the Higgs and top are boosted works well to distinguish the
signal and backgrounds. Such region can be characterized by the angular
distances $\Delta R\leq 1.4$ of the leptons from the Higgs boson and $\Delta
R>1.8$ between the final jet from top quark and leptons from the Higgs boson.
And these cuts work better on 100TeV hadron collider. We also checked the
possibility to search this process at the hadron colliders. We found that the
$WW^{*}$ channel is a promising channel on both 14TeV and 100TeV colliders
while the cross-section of $ZZ^{*}$ channel seems a little small. We combined
different channels and obtained the expected upper limits of the couplings.
These limits can be presented as $Br(t\to Hc)\leq 4.09\times 10^{-3}$ and
$Br(t\to Hu)\leq 6.06\times 10^{-4}$. We can see that the process $pp\to
t/\bar{t}H$ is not as good as $t\to Hq$ in setting the limits of $y_{tq}$.
However, this process can be used to distinguish the $y_{tq}$ couplings, if
the FCNC decay of the top quark $t\to Hq$ is observed.

###### Acknowledgements.

This work is supported in part by the Natural Science Foundation of Shandong
Province under grant No. ZR2020MA094.

## References

  * [1] G. Aad et al. [ATLAS], Phys. Lett. B 716, 1-29 (2012) doi:10.1016/j.physletb.2012.08.020 [arXiv:1207.7214 [hep-ex]].
  * [2] S. Chatrchyan et al. [CMS], Phys. Lett. B 716, 30-61 (2012) doi:10.1016/j.physletb.2012.08.021 [arXiv:1207.7235 [hep-ex]].
  * [3] B. Mele, S. Petrarca and A. Soddu, Phys. Lett. B 435, 401-406 (1998) doi:10.1016/S0370-2693(98)00822-3 [arXiv:hep-ph/9805498 [hep-ph]].
  * [4] G. Eilam, J. L. Hewett and A. Soni, Phys. Rev. D 44, 1473-1484 (1991) [erratum: Phys. Rev. D 59, 039901 (1999)] doi:10.1103/PhysRevD.44.1473
  * [5] Y. L. Wu and L. Wolfenstein, Phys. Rev. Lett. 73, 1762-1764 (1994) doi:10.1103/PhysRevLett.73.1762 [arXiv:hep-ph/9409421 [hep-ph]].
  * [6] D. Atwood, L. Reina and A. Soni, Phys. Rev. D 55, 3156-3176 (1997) doi:10.1103/PhysRevD.55.3156 [arXiv:hep-ph/9609279 [hep-ph]].
  * [7] S. S. Bao, F. Su, Y. L. Wu and C. Zhuang, Phys. Rev. D 77, 095004 (2008) doi:10.1103/PhysRevD.77.095004 [arXiv:0801.2596 [hep-ph]].
  * [8] H. E. Haber and G. L. Kane, Phys. Rept. 117, 75-263 (1985) doi:10.1016/0370-1573(85)90051-1
  * [9] J. Cao, C. Han, L. Wu, J. M. Yang and M. Zhang, Eur. Phys. J. C 74, no.9, 3058 (2014) doi:10.1140/epjc/s10052-014-3058-1 [arXiv:1404.1241 [hep-ph]].
  * [10] A. Dedes, M. Paraskevas, J. Rosiek, K. Suxho and K. Tamvakis, JHEP 11, 137 (2014) doi:10.1007/JHEP11(2014)137 [arXiv:1409.6546 [hep-ph]].
  * [11] J. A. Aguilar-Saavedra, Acta Phys. Polon. B 35, 2695-2710 (2004) [arXiv:hep-ph/0409342 [hep-ph]].
  * [12] S. S. Bao, Commun. Theor. Phys. 71, no.9, 1093-1096 (2019) doi:10.1088/0253-6102/71/9/1093
  * [13] J. B. Guimarães da Costa et al. [CEPC Study Group], [arXiv:1811.10545 [hep-ex]].
  * [14] K. F. Chen, W. S. Hou, C. Kao and M. Kohda, Phys. Lett. B 725, 378-381 (2013) doi:10.1016/j.physletb.2013.07.060 [arXiv:1304.8037 [hep-ph]].
  * [15] C. Kao, H. Y. Cheng, W. S. Hou and J. Sayre, Phys. Lett. B 716, 225-230 (2012) doi:10.1016/j.physletb.2012.08.032 [arXiv:1112.1707 [hep-ph]].
  * [16] D. Atwood, S. K. Gupta and A. Soni, JHEP 10, 057 (2014) doi:10.1007/JHEP10(2014)057 [arXiv:1305.2427 [hep-ph]].
  * [17] T. Han, J. Jiang and M. Sher, Phys. Lett. B 516, 337-344 (2001) doi:10.1016/S0370-2693(01)00949-2 [arXiv:hep-ph/0106277 [hep-ph]].
  * [18] L. Wu, JHEP 02, 061 (2015) doi:10.1007/JHEP02(2015)061 [arXiv:1407.6113 [hep-ph]].
  * [19] Y. Wang, F. P. Huang, C. S. Li, B. H. Li, D. Y. Shao and J. Wang, Phys. Rev. D 86, 094014 (2012) doi:10.1103/PhysRevD.86.094014 [arXiv:1208.2902 [hep-ph]].
  * [20] A. Greljo, J. F. Kamenik and J. Kopp, JHEP 07, 046 (2014) doi:10.1007/JHEP07(2014)046 [arXiv:1404.1278 [hep-ph]].
  * [21] D. Atwood, S. K. Gupta and A. Soni, JHEP 04, 035 (2013) doi:10.1007/JHEP04(2013)035 [arXiv:1301.2250 [hep-ph]].
  * [22] O. M. Ozsimsek, V. Ari and O. Cakir, [arXiv:2204.13139 [hep-ph]].
  * [23] [ATLAS], ATLAS-CONF-2022-014.
  * [24] A. Tumasyan et al. [CMS], JHEP 02, 169 (2022) doi:10.1007/JHEP02(2022)169 [arXiv:2112.09734 [hep-ex]].
  * [25] J. Alwall, R. Frederix, S. Frixione, V. Hirschi, F. Maltoni, O. Mattelaer, H. S. Shao, T. Stelzer, P. Torrielli and M. Zaro, JHEP 07, 079 (2014) doi:10.1007/JHEP07(2014)079 [arXiv:1405.0301 [hep-ph]].
  * [26] C. Bierlich, S. Chakraborty, N. Desai, L. Gellersen, I. Helenius, P. Ilten, L. Lönnblad, S. Mrenna, S. Prestel and C. T. Preuss, et al. [arXiv:2203.11601 [hep-ph]].
  * [27] J. de Favereau et al. [DELPHES 3], JHEP 02, 057 (2014) doi:10.1007/JHEP02(2014)057 [arXiv:1307.6346 [hep-ex]].
  * [28] R. D. Ball, V. Bertone, S. Carrazza, C. S. Deans, L. Del Debbio, S. Forte, A. Guffanti, N. P. Hartland, J. I. Latorre and J. Rojo, et al. Nucl. Phys. B 867, 244-289 (2013) doi:10.1016/j.nuclphysb.2012.10.003 [arXiv:1207.1303 [hep-ph]].
  * [29] J. M. Butterworth, I. Ochoa and T. Scanlon, Eur. Phys. J. C 75, no.8, 366 (2015) doi:10.1140/epjc/s10052-015-3592-5 [arXiv:1506.04973 [hep-ph]].
  * [30] M. Dittmar and H. K. Dreiner, Phys. Rev. D 55, 167-172 (1997) doi:10.1103/PhysRevD.55.167 [arXiv:hep-ph/9608317 [hep-ph]].
  * [31] G. Cowan, K. Cranmer, E. Gross and O. Vitells, Eur. Phys. J. C 71, 1554 (2011) [erratum: Eur. Phys. J. C 73, 2501 (2013)] doi:10.1140/epjc/s10052-011-1554-0 [arXiv:1007.1727 [physics.data-an]].
  * [32] [ATLAS], ATL-PHYS-PUB-2016-019.


# Inertial active harmonic particle with memory escape induced by viscoelastic
suspension

F Adersh  Department of Physics, University of Kerala, Kariavattom,
Thiruvananthapuram-$695581$, India  M Muhsin  Department of Physics,
University of Kerala, Kariavattom, Thiruvananthapuram-$695581$, India  M Sahoo
<EMAIL_ADDRESS>Department of Physics, University of Kerala, Kariavattom,
Thiruvananthapuram-$695581$, India

###### Abstract

We investigate the self-propulsion of an inertial active particle confined in
a two-dimensional harmonic trap. The particle is suspended in a non-Newtonian
or viscoelastic suspension with a friction kernel that decays exponentially
with a time constant characterizing the memory timescale or transient
elasticity of the medium. By solving the associated non-Markovian dynamics, we
identify two regimes in parameter space distinguishing the oscillatory and
non-oscillatory behavior of the particle motion. By simulating the particle
trajectories and exactly calculating the steady state probability distribution
functions and mean square displacement, interestingly, we observe that with an
increase in the memory time scale, the elastic bound of suspension dominates
over the influence of harmonic trap. As a consequence, the particle can escape
out of the trap without approaching steady state. On the other hand, with an
increase in the duration of the activity, the particle becomes trapped by the
harmonic confinement.

##  I INTRODUCTION

The physics of active matter represents a rapidly advancing field of research
that has gained substantial interest across various scientific disciplines [1,
2, 3, 4, 5], particularly in the realm of biophysics and bioengineering.
Active or self-propelling systems constitute entities that are driven out of
equilibrium by harnessing energy from their environment to generate directed
motion. Active matter primarily includes microorganisms like E.coli, ciliates,
and other motile bacteria, but they can also be synthesized artificially.
Janus particles [6, 7] and microrobots [8] are some of the prominent examples
of artificially synthesized active matter. They can mimic biological motility
and find applications in many emerging scientific and technological domains
like material science and nanotechnology. Some of the notable theoretical
models for exploring active matter include active Ornstein-Uhlenbeck particle
(AOUP) model[9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20], Active Brownian
particle model (ABP) [21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34],
run and tumble particle (RTP) model [35, 36], etc. The AOUP model introduces
activity into the system through a stochastic force that follows the Ornstein-
Uhlenbeck process. The ABP model considers both translational and rotational
diffusion of particles, while RTPs exhibit a “run” phase characterized by
ballistic motion with constant speed followed by periodic tumbling.

An inertial active particle while self-propelling within a harmonic
confinement can exhibit oscillatory motion and the frequency of oscillation
depends on the inertia and strength of the harmonic confinement [37]. Such an
active harmonic particle while self-propelling in a two-dimensional (2d)
plane, the oscillatory motion turns in the performance of rotational
trajectories [13, 38, 39]. Further, an active particle can also exhibit
oscillation in its motion while suspended in a viscoelastic environment [40,
41]. Such induced oscillations often serve as distinctive characteristic of
non-equilibrium particle activity [42]. Moreover, such a viscoelastic particle
while confined in a plane, it performs circular motion [43].

In this paper, we delve into the dynamics of an inertial active Ornstein-
Uhlenbeck particle (AOUP) suspended in a non-Newtonian bath and confined
within a harmonic potential. The viscoelastic memory of the environment
influences the particle motion in both qualitative and quantitative ways.
Examples include the memory induced delay between the effective self-
propulsion force and particle orientation in an overdamped active particle
[44], strong enhancement of directional changes in the periodically modulated
Run and Tumble particles in viscoelastic bath [45], effective repulsion among
active particles induced by viscoelastic fluid while moving near a
boundary[46], etc. By simulating and exactly solving our present model
dynamics, we observe that for specific choice of parameters, the particle
exhibits oscillatory motion and the frequency of oscillation exhibits a non-
monotonic dependence on both inertial and harmonic time scales. The simulated
particle trajectories confirm that the particle adopts rotational motion in
the oscillatory regime of parameter space. Moreover, an increase in the memory
timescale of the medium results in the effective dominance of elastic
dissipation over harmonic confinement. As a consequence, the particle can
escape out of the potential or blow off without approaching a steady state. On
the other hand, an increase in the activity time results in trapping of the
particle by the harmonic confinement. These findings are consistent with the
analytically computed steady state results for probability distribution
functions and mean square displacement.

##  II MODEL

We consider the 2d motion of an inertial active Ornstein-Uhlenbeck particle in
a viscoelastic environment. The particle is confined in a harmonic potential
$U(x,y)=\frac{1}{2}k(x^{2}+y^{2})$ with $k$ as the harmonic constant. The
dynamics is non-Markovian due to the non-Newtonian nature of the surrounding
medium and hence can be described using the generalized Langevin’s equation of
motion [41, 20]

$m\ddot{\bf r}(t)=-\gamma\int_{0}^{t}\lambda(t-t^{\prime})\dot{\bf
r}(t^{\prime})\,dt^{\prime}-\boldsymbol{\nabla}U+\boldsymbol{\xi}(t).$ (1)

Here, ${\bf r}(t)=x(t)\hat{i}+y(t)\hat{j}$ is the position vector of the
particle in the $x-y$ plane and $m$ represents the mass of the particle. The
first term in RHS defines the viscous drag and is characterized by an
exponentially decaying memory kernel $\lambda(\tau)$ of the form

$\lambda\left(\tau\right)=\begin{cases}\frac{1}{t_{c}^{\prime}}e^{-\frac{\tau}{t_{c}^{\prime}}};&\tau\geq
0,\\\ 0~{}~{}~{}~{}~{}~{}~{}~{}~{};&\tau<0.\end{cases}$ (2)

The above memory kernel represents that the medium exhibits a transient
elasticity that exponentially decays to viscous behavior with a timescale
$t_{c}^{\prime}$. Hence, $t_{c}^{\prime}$ can be regarded as the memory time
scale or elastic dissipation timescale of the system. The term
$\boldsymbol{\xi}(t)$ in Eq. (1) represents the athermal noise associated with
the dynamics which follows the Ornstein-Uhlenbeck (OU) process

$t_{c}\dot{\boldsymbol{\xi}}(t)=-\boldsymbol{\xi}(t)+\sqrt{D}\boldsymbol{\eta}(t).$
(3)

Here, $t_{c}$ represents the self-propulsion or activity timescale of the
dynamics, $D$ being the strength of noise and $\boldsymbol{\eta}(t)$ is the
delta-correlated white noise. $\boldsymbol{\xi}(t)$ has the properties

$\langle\xi_{\alpha}(t)\rangle=0,\qquad\langle\xi_{\alpha}(t)\xi_{\beta}(t^{\prime})\rangle=\frac{\delta_{\alpha\beta}}{2t_{c}}e^{-\frac{|t-t^{\prime}|}{t_{c}}}.$
(4)

When $t_{c}=t_{c}^{\prime}$, the generalized fluctuation-dissipation theorem
is validated and the system approaches thermal equilibrium for $D=2\gamma
k_{B}T$.

Now we define the physical quantities of interest. The mean displacement (MD)
of the particle from the initial position is given by

$\langle\Delta{\bf r}(t)\rangle=\langle{\bf r}(t)-{\bf r}(0)\rangle.$ (5)

Similarly, the mean square displacement (MSD) can be obtained as

$\begin{split}\langle{\Delta{\bf r}(t)}^{2}\rangle&=\langle\left({\bf
r}(t)-{\bf r}(0)\right)^{2}\rangle.\end{split}$ (6)

To obtain the particle trajectory, we perform a numerical simulation of
Langevin dynamics [Eq. (1)]. The integration of the equation of motion [Eq.
(1)] is carried out using a second-order modified Euler method with a timestep
of $10^{-3}$. The OU noise is realized using the Fox algorithm[47]. In the
next section, we discuss the results obtained from the numerical simulation
and analytical calculations.

##  III RESULTS AND DISCUSSION

Figure 1: (a) Phase diagram separating the oscillatory and non-oscillatory
regimes in the $t_{c}^{\prime}\Gamma-t_{c}^{\prime}\omega_{0}$ parameter space
[Eq. (7)]. The color map in (a) shows the evolution of frequency of
oscillation ($\nu$) with $t_{c}^{\prime}\Gamma$ and
$t_{c}^{\prime}\omega_{0}$. (b) $\nu$ vs $\Gamma$ for different values of
$\omega_{0}$ and for a fixed $t_{c}^{\prime}=1$. (c) $\nu$ vs $\omega_{0}$ for
different values of $\Gamma$ and for a fixed $t_{c}^{\prime}=1$. (d) The time
evolution of the x-component of mean displacement (MD) $\langle\Delta
x(t)\rangle$ for different values of $\omega_{0}$ and $\Gamma$ with
$t_{c}^{\prime}=1.0$. The other common parameters are $m=t_{c}=D=1$.

.

Solving the dynamics [Eq. (1)], with initial conditions ${\bf r}(0)={\bf
r}_{0}$ and ${\bf\dot{r}}(0)={\bf v}_{0}$, the solution can be obtained as

${\bf r(t)}=\sum_{i=1}^{3}a_{i}\biggr{[}e^{s_{i}t}\left({\bf v_{0}}+s_{i}{\bf
r_{0}}\right)+\int_{0}^{t}e^{s_{i}(t-t^{\prime})}\xi(t^{\prime})\,dt^{\prime}\biggr{]}.$
(7)

Here, $s_{i}$’s are the roots of the equation

$s^{3}+\frac{s^{2}}{t_{c}^{\prime}}+s\left(\frac{\Gamma}{t_{c}^{\prime}}+\omega_{0}^{2}\right)+\frac{\omega_{0}^{2}}{t_{c}^{\prime}}=0,$
(8)

which are given by

$s_{i}=\frac{1}{3t_{c}^{\prime}}\left[-1-\frac{\omega^{i-1}\Delta_{0}}{\sqrt[3]{\Delta_{1}}}+\omega^{i-1}\sqrt[3]{\Delta_{1}}\right],\quad
i\in\\{1,2,3\\}$ (9)

with

$\begin{split}\Delta_{0}&=-1+3\left[t_{c}^{\prime}\Gamma+(t_{c}^{\prime}\omega_{0})^{2}\right],\\\
\Delta_{1}&=\frac{-2+9t_{c}^{\prime}\Gamma-18(t_{c}^{\prime}\omega_{0})^{2}+\sqrt{-\Delta}}{2},\\\
\Delta&=-\left[2-9t_{c}^{\prime}\Gamma+18(t_{c}^{\prime}\omega_{0})^{2}\right]^{2}-4\Delta_{0}^{3},\quad\text{and}\\\
\omega&=\frac{-1+j\sqrt{3}}{2}.\end{split}$ (10)

Here, $\Gamma=\frac{\gamma}{m}$, $\omega_{0}=\sqrt{\frac{k}{m}}$ and
$j=\sqrt{-1}$. The coefficients $a_{i}$’s in Eq. (7) are given by

$\begin{split}&a_{1}=\frac{s_{1}t_{c}^{\prime}+1}{\left(s_{1}-s_{2}\right)\left(s_{1}-s_{3}\right)t_{c}^{\prime}},~{}a_{2}=\frac{s_{2}t_{c}^{\prime}+1}{\left(s_{2}-s_{1}\right)\left(s_{2}-s_{3}\right)t_{c}^{\prime}}\\\
&\text{and}~{}~{}a_{3}=\frac{s_{3}t_{c}^{\prime}+1}{\left(s_{2}-s_{3}\right)\left(s_{1}-s_{3}\right)t_{c}^{\prime}}.\end{split}$
(11)

Solving Eq. (8), it is realized that the roots of Eq. (8) can be real as well
as complex. Positive and negative values of $\Delta$ result in some of the
roots of Eq. (8) to be complex. The complex roots are responsible for the
oscillatory solution of the dynamics [Eq. (7)], whereas the real roots provide
an exponentially decaying solution without any oscillatory behavior. Thus,
there exist two separate regimes associated with the oscillatory and non-
oscillatory behavior of the solution. Since the roots $s_{i}$ are independent
of $t_{c}$, these oscillatory and non-oscillatory regimes are presented in a
phase diagram of $t_{c}^{\prime}\omega_{0}$-$t_{c}^{\prime}\Gamma$ parameter
space [see Fig. 1(a)]. The boundary separating these two regimes in
$t_{c}^{\prime}\omega_{0}$ – $t_{c}^{\prime}\Gamma$ space is given by
$\Delta=0$. The color map in Fig. 1(a) shows the evolution of frequency of
oscillation ($\nu$) with both $\omega_{0}$ and $\Gamma$ for a fixed
$t_{c}^{\prime}$ value. The frequency of oscillation is computed as the
imaginary part of $s_{i}$ (Eq. (9)). In Fig. 1(b), we have shown the variation
of $\nu$ with $\Gamma$ for different values of $\omega_{0}$ by keeping
$t_{c}^{\prime}$ fixed. For small values of $\omega_{0}$, $\nu$ as a function
of $\Gamma$ initially decreases, becomes zero and shows a plateau for the
range of $\Gamma$ that falls in the non-oscillatory regime of Fig. 1(a) and
finally increases with increase in $\Gamma$ value. For higher values of
$\omega_{0}$, this non-monotonic behavior of $\nu$ with $\Gamma$ disappears
and it monotonically increases with $\Gamma$. Similarly, Fig. 1(c) shows the
variation of $\nu$ with $\omega_{0}$ for different values of $\Gamma$ and for
a fixed $t_{c}^{\prime}$. For lower $\Gamma$ values, $\nu$ as a function of
$\omega_{0}$ shows a non-monotonic behavior, it decreases, approaches zero,
and finally increases with $\omega_{0}$. For higher $\Gamma$ values, $\nu$
becomes an increasing function of $\omega_{0}$. For further exploration of
these two regimes, in Fig. 1(d), we plot the x-component of MD $\langle\Delta
x(t)\rangle$ for different values of $\Gamma$ and $\omega_{0}$ keeping
$t_{c}^{\prime}=1$. Using Eq. (7), MD [Eq. (5)] of the particle can be
calculated as

$\langle\Delta{\bf r}(t)\rangle=\sum_{i=1}^{3}a_{i}e^{s_{i}t}\Bigr{[}{\bf
v_{0}}+{\bf r_{0}}s_{i}\Bigr{]}-{\bf r_{0}}.$ (12)

In the lower powers of $t$, it can be expanded as

$\langle\Delta{\bf r}(t)\rangle=-{\bf r_{0}}+{\bf v_{0}}t-\frac{{\bf
r_{0}}\left(\Gamma+4t_{c}^{\prime}\omega_{0}^{2}\right)}{8t_{c}^{\prime}}t^{2}+O\left(t^{3}\right).$
(13)

From Fig. 1(d), it is to be noted that $\langle\Delta x(t)\rangle$ does not
show any oscillatory behavior for the parameters that lie in the no-
oscillation regimes of $t_{c}^{\prime}\Gamma$-$t_{c}^{\prime}\omega_{0}$
parameter space (solid blue curve for $\Gamma=0.2$ and $\omega_{0}=0.1$). On
the other hand, $\langle\Delta x(t)\rangle$ exhibits an intermediate time
oscillatory behavior(curves other than blue in color) for the parameters that
fall in the oscillatory regime of
$t_{c}^{\prime}\Gamma$-$t_{c}^{\prime}\omega_{0}$ parameter space.

Figure 2: The simulated particle trajectories in the $xy$-plane for different
values of $\omega_{0}$ are plotted for $t_{c}=1$ in (a),(b) and (c), for
$t_{c}=5$ in (d), (e) and (f), and for $t_{c}=10$ in (g), (h) and (i),
respectively. The other common parameters are $m=t_{c}^{\prime}=D=1$, and
$\Gamma=0.2$. The color map shows the time evolution of the trajectories.

In Fig. 2, we plot the instantaneous 2d particle trajectories for different
values of $\omega_{0}$ and $t_{c}$. Each column of Fig. 2 corresponds to a
fixed value $\omega_{0}$ and each row corresponds to a fixed value of $t_{c}$.
For a given set of parameters, with an increase in $\omega_{0}$ value from
left to right in any of the rows, the particle passes from no-oscillatory to
oscillatory regime. The trajectories reflect random self-propulsion for the
parameters that lie in the no-oscillatory regimes [see Figs. 2(a), (d) and
(g)]. However, the particle performs rotational trajectories for the
parameters that lie in the oscillatory regime of parameter space [see Figs.
2(b), (c), (e), (f), (h) and (i)]. Further, the trajectories show stronger
confinement of the particle with increase in $\omega_{0}$ value as expected.
Similarly, in any of the columns, the trajectories get suppressed with
increase in $t_{c}$ value. This observation suggests the trapping or
confinement of the particle around the centre of the potential with increase
in persistent duration of activity.

Figure 3: The simulated particle trajectories in the $xy$-plane for different
values of $\omega_{0}$ are plotted for $t_{c}^{\prime}=0.1$ in (a), (b) and
(c), for $t_{c}^{\prime}=1$ in (d), (e) and (f), and for $t_{c}^{\prime}=5.0$
in (g), (h) and (i), respectively. The other common parameters are $m=1$,
$\Gamma=0.05$, $t_{c}=1.0$ and $D=1.0$. The color map shows the time evolution
of the trajectories.

Similarly, the 2d particle trajectories for different values of $\omega_{0}$
and $t_{c}^{\prime}$ are shown in Fig. 3, with $\omega_{0}$ increasing along a
row from left to right and $t_{c}^{\prime}$ increasing along a column from top
to bottom. Similar to the observations in Fig. 2, for a fixed $t_{c}^{\prime}$
value, with increase in $\omega_{0}$ along a row, the parameters are chosen
such that the particle passes from no-oscillatory to oscillatory regime. The
particle performs random self-propulsion in the no-oscillatory regime, i.e.,
for $\Delta>0$ [see Fig. 3(a), (d) and (g)]. For higher $\omega_{0}$ values,
the particle enters the oscillatory regime (since $\Delta$ becomes negative)
and makes a transition from its random activity to twisting or rotational
motion [see Figs. 3(b), (c), (e), (f), (h) and (i)]. Further, with increase in
$t_{c}^{\prime}$ value from top to bottom in a column, the trajectories get
enhanced. At the same time, the trajectories become more asymmetric and take
almost elliptical shape. The enhancement of trajectory as an increasing
function of $t_{c}^{\prime}$ can be attributed to the fact that for high value
of $t_{c}^{\prime}$, the system retains memory for a prolonged period,
resulting in less frequent change in the direction of the particle velocity.
As a consequence, the particle takes longer time to complete the orbit.

Next, we have calculated the MSD [Eq. (6)] using Eq. (7) as

$\begin{split}\langle\Delta{\bf r}^{2}(t)\rangle=&\langle\left({\bf r}(t)-{\bf
r}(0)\right)^{2}\rangle=\left(\langle\Delta{\bf
r}(t)\rangle\right)^{2}+\sum_{i=1}^{3}\sum_{j=1}^{3}\frac{a_{i}a_{j}^{*}D}{m^{2}}\Biggl{[}\frac{t_{c}e^{t\left(s^{*}_{j}-\frac{1}{t_{c}}\right)}}{\left(t_{c}s_{i}+1\right)\left(1-t_{c}s^{*}_{j}\right)}+\frac{t_{c}e^{t\left(s_{i}-\frac{1}{t_{c}}\right)}}{\left(1-t_{c}s_{i}\right)\left(t_{c}s^{*}_{j}+1\right)}\\\
&+\frac{t_{c}\left(s_{i}+s^{*}_{j}\right)-2}{\left(s_{i}+s^{*}_{j}\right)\left(t_{c}s_{i}-1\right)\left(t_{c}s^{*}_{j}-1\right)}+\frac{e^{t\left(s_{i}+s^{*}_{j}\right)}\left(t_{c}\left(s_{i}+s^{*}_{j}\right)+2\right)}{\left(s_{i}+s^{*}_{j}\right)\left(t_{c}s_{i}+1\right)\left(t_{c}s^{*}_{j}+1\right)}\Biggr{]}.\end{split}$
(14)

Figure 4: The MSD [Eq. (14)] as a function of $t$ for different values of
$\Gamma$ in (a) for $t_{c}=t_{c}^{\prime}=1.0$ and $\omega_{0}=0.1$, for
different values of $\omega_{0}$ in (b) for $t_{c}=t_{c}^{\prime}=1.0$ and
$\Gamma=0.2$, for different values of $t_{c}$ in (c) for a fixed
$t_{c}^{\prime}=1$, and for different values of $t_{c}^{\prime}$ in (d) for a
fixed $t_{c}=1$, respectively. The insets of (a) and (b) show the exponent
$\alpha$ as a function of $t$ corresponding to the respective MSD plots. The
dashed horizontal lines in (c) and (d) correspond to the steady state
equilibrium value of MSD. The other common parameters are $D=m=1.$

.

Figure 5: The 2d plot of steady state MSD $\langle\Delta r^{2}\rangle_{s}$
[Eq. (28)] as a function of both $t_{c}$ and $t_{c}^{\prime}$. The other
common parameters are $\Gamma=1.0$, $m=1$, $D=1$ and $\omega_{0}=1.0$.  Figure
6: The steady state position probability distribution [Eq. (32)] for different
values of $t_{c}^{\prime}$ fixing $t_{c}(=5.0)$ are shown in (a)-(c) and for
different values of $t_{c}$ for a fixed $t_{c}^{\prime}(=5)$ is shown in
(d)-(f). Other parameters are $m=\Gamma=\omega_{0}=D=1$.

In Fig. 4, we present $\langle\Delta{\bf r}^{2}(t)\rangle$ as a function of
$t$ for different values of $\Gamma$ in Fig. 4(a), $\omega_{0}$ in Fig. 4 (b),
$t_{c}$ in Fig. 4(c), and $t_{c}^{\prime}$ in Fig. 4(d), respectively. In
order to explore different time regimes of motion, we introduce a parameter
$\alpha$ such that $\langle\Delta{\bf r}^{2}(t)\rangle\propto t^{\alpha}$,
i.e.,

$\alpha=\frac{\partial\log\langle\Delta{\bf r}^{2}(t)\rangle}{\partial\log
t}.$ (15)

The insets of Fig. 4(a) and (b) show the variation of the exponent $\alpha$
with $t$. The initial transient time MSD is always found to be ballistic
(since $\langle\Delta{\bf r}^{2}(t)\rangle\propto t^{2}$ or $\alpha=2$) and
the steady state or long time regime is always non-diffusive (since
$\langle\Delta{\bf r}^{2}(t)\rangle$ is independent of $t$ or $\alpha=0$).
However, the intermediate time regimes of MSD is found to be oscillatory for
the parameters that fall in the oscillatory regime of parameter space. In Fig.
4(a), the variation of MSD as a function of $t$ for low and high values of
$\Gamma$ ($\Gamma=0.05,0.5,1.0$) show intermediate time oscillations and that
is why the exponent $\alpha$ as a function of $t$ shows intermediate time
oscillations for $\Gamma=0.05,0.5$, and $1.0$, respectively, whereas the MSD
is non-oscillatory for $\Gamma=0.2$. With increase in $\omega_{0}$ value in
Fig. 4 (b), the system makes a transitions from no-oscillatory phase to
oscillatory phase and results intermediate time oscillations in MSD for larger
values of $\omega_{0}$ (for $\omega_{0}=0.5$ and $\omega_{0}=1$). The same can
also be confirmed from the inset of Fig. 4 (b) which shows the intermediate
time oscillatory behaviour of $\alpha$ for $\omega_{0}=0.5$ and
$\omega_{0}=1$. Figures 4 (c) and (d) show the variation of MSD as a function
of $t$ for different values of $t_{c}$ and $t_{c}^{\prime}$, respectively.
With an increase in $t_{c}$ values, the steady state MSD decreases (see Fig. 4
(c)) and it increases with an increase in $t_{c}^{\prime}$ values (see Fig. 4
(d)). At the same time, the initial ballistic regime of MSD gets reduced with
increase in $t_{c}$ value and increases with increase in $t_{c}^{\prime}$
values.

We exactly evaluate the steady state MSD and probability distribution function
using steady state correlation matrix method formalism [13, 48]. Because of
the linear dynamics of the model and Gaussian nature of noise, the steady
state probability distribution is Gaussian. In order for evaluating these two
quantities, we introduce the correlation matrix $\Xi$ with components
$\chi_{i}$, $\chi_{j}$ such that

$[\Xi_{i,j}]=\langle\chi_{i}\chi_{j}\rangle-\langle\chi_{i}\rangle\langle\chi_{j}\rangle$
(16)

and the distribution function $\mathcal{P}({\boldsymbol{\chi}})$ as

$\mathcal{P}({\boldsymbol{\chi}})=\frac{1}{(2\pi)^{4}\sqrt{Det(\Xi)}}e^{\left[{-\frac{1}{2}\left(\boldsymbol{\chi}^{T}\right)\Xi^{-1}\left(\boldsymbol{\chi}\right)}\right]}.$
(17)

Here, $\boldsymbol{\chi}$ is the column vector whose elements
$\chi_{i}\in\left(x\ y\ v_{x}\ v_{y}\ w_{x}\ w_{y}\ \xi_{x}\ \xi_{y}\right)$
with $v_{x}$ and $v_{y}$ being the $x$ and $y$ components of velocity. The
variables $w_{x}$ and $w_{y}$ are the $x$ and $y$ components of the vector
${\bf w}$ given by

${\bf
w}=\Gamma\int_{0}^{t}\lambda\left(t-t^{\prime}\right){\bf\dot{r}}(t^{\prime})\,dt^{\prime}.$
(18)

Introduction of ${\bf w}$ splits the non-Markovian model dynamics [Eq. (1)]
into a set of Markovian equations as follows

$\displaystyle{\bf\dot{r}}$ $\displaystyle={\bf v}$ (19)
$\displaystyle{\bf\dot{v}}$ $\displaystyle=-{\bf w}-\omega_{0}^{2}{\bf
r}+\frac{\sqrt{D}}{m}\boldsymbol{\xi}$ (20) $\displaystyle{\bf\dot{w}}$
$\displaystyle=-\frac{{\bf
w}}{t_{c}^{\prime}}+\frac{\Gamma}{t_{c}^{\prime}}{\bf v}$ (21)
$\displaystyle{\bf\dot{\boldsymbol{\xi}}}$
$\displaystyle=-\frac{\boldsymbol{\xi}}{t_{c}}+\frac{\boldsymbol{\eta}}{t_{c}}.$
(22)

The above set of equations can be expressed as

$\dot{\boldsymbol{\chi}}=A\boldsymbol{\chi}+B\boldsymbol{\eta^{\prime}}$ (23)

Here, the matrices $A$, $B$ and $\boldsymbol{\eta^{\prime}}$ are given by

$A=\begin{pmatrix}0&0&1&0&0&0&0&0\\\ 0&0&0&1&0&0&0&0\\\
-\omega_{0}^{2}&0&0&0&-1&0&\frac{\sqrt{D}}{m}&0\\\
0&-\omega_{0}^{2}&0&0&0&-1&0&\frac{\sqrt{D}}{m}\\\
0&0&\frac{\Gamma}{t_{c}^{\prime}}&0&\frac{-1}{t_{c}^{\prime}}&0&0&0\\\
0&0&0&\frac{\Gamma}{t_{c}^{\prime}}&0&\frac{-1}{t_{c}^{\prime}}&0&0\\\
0&0&0&0&0&0&\frac{-1}{t_{c}}&0\\\
0&0&0&0&0&0&0&\frac{-1}{t_{c}}\end{pmatrix},$ (24)
$B=\begin{pmatrix}0&0&0&0&0&0&0&0\\\ 0&0&0&0&0&0&0&0\\\ 0&0&0&0&0&0&0&0\\\
0&0&0&0&0&0&0&0\\\ 0&0&0&0&0&0&0&0\\\ 0&0&0&0&0&0&0&0\\\
0&0&0&0&0&0&\frac{1}{t_{c}}&0\\\ 0&0&0&0&0&0&0&\frac{1}{t_{c}}\\\
\end{pmatrix},$ (25)

and

$\boldsymbol{\eta^{\prime}}=(0\ 0\ 0\ 0\ 0\ 0\ \eta_{x}\ \eta_{y})^{T}.$ (26)

As per the correlation matrix formalism, the correlation matrix $\Xi$ can be
shown to satisfy the equation

$A\cdot\Xi+\Xi\cdot A^{T}+B\ B^{T}=0.$ (27)

Finally, solving the above Eq. (27), the MSD at steady state
$\langle\Delta{\bf r}^{2}\rangle_{s}$ can be calculated as

$\begin{split}\langle\Delta{\bf r}^{2}\rangle_{s}&=\Xi_{1,1}+\Xi_{2,2}\\\
&=\frac{D\left[t_{c}+t_{c}^{\prime}+t_{c}^{2}\Gamma+t_{c}^{\prime
2}(t_{c}+t_{c}^{\prime})\omega_{0}^{2}\right]}{m^{2}\Gamma\omega_{0}^{2}\left[t_{c}^{\prime}+t_{c}(1+t_{c}\Gamma+t_{c}(t_{c}+t_{c}^{\prime})\omega_{0}^{2})\right]}\\\
&=\frac{D}{\Gamma
m^{2}\omega_{0}^{2}}+\frac{D\left(t_{c}^{\prime}-t_{c}\right)\left(t_{c}^{\prime}+t_{c}\right){}^{2}}{\Gamma
m^{2}\left(\Gamma
t_{c}^{2}+\omega_{0}^{2}t_{c}^{2}t_{c}^{\prime}+\omega_{0}^{2}t_{c}^{3}+t_{c}^{\prime}+t_{c}\right)}.\end{split}$
(28)

From the above equation, it is confirmed that for $t_{c}=t_{c}^{\prime}$
limit, $\langle\Delta{\bf r}^{2}\rangle_{s}$ reduces to $\frac{D}{\Gamma
m^{2}\omega_{0}^{2}}$. For $D=2\gamma k_{B}T$, it approaches the equilibrium
value given by $\langle\Delta{\bf r}^{2}\rangle_{eq}=\frac{2k_{B}T}{k}$.
Similarly, in $\omega_{0}\to\infty$ limit, $\langle\Delta{\bf
r}^{2}\rangle_{s}$ approaches the equilibrium value for $D=2\gamma k_{B}T$
even when $t_{c}\neq t_{c}^{\prime}$. However, in this limit, the system
doesn’t approach equilibrium. The system approaches thermal equilibrium for
$D=2\gamma k_{B}T$ in $t_{c}=t_{c}^{\prime}$ limit and the equilibrium value
follows the equipartition theorem

$\frac{1}{2}k\langle\Delta{\bf r}^{2}\rangle_{eq}=k_{B}T.$ (29)

The first term of Eq. (28) is independent of $t_{c}$ and $t_{c}^{\prime}$.
With increase in $t_{c}^{\prime}$ value, the second term of Eq. (28) becomes
dominant and hence in $t_{c}^{\prime}\to\infty$ limit, this term diverges. As
a result, the steady state MSD diverges, i.e.,

$\lim\limits_{t_{c}^{\prime}\to\infty}\langle\Delta{\bf
r}^{2}\rangle_{s}=\infty.$ (30)

Similarly, with increase in $t_{c}$ value, the second term of Eq. (28) becomes
negative, as a result the $\langle\Delta{\bf r}^{2}\rangle_{s}$ decreases. In
$t_{c}\to\infty$ limit, the second term of Eq. (28) becomes $-\frac{D}{\Gamma
m^{2}\omega_{0}^{2}}$, which is exactly equal to the negative of the first
term. Hence, in this limit, the steady state MSD vanishes, i.e.,

$\lim\limits_{t_{c}\to\infty}\langle\Delta{\bf r}^{2}\rangle_{s}=0.$ (31)

These results are summarized in Fig. 5, where we present the 2d plot of
$\langle\Delta{\bf r}^{2}\rangle_{s}$ as a function of both $t_{c}$ and
$t_{c}^{\prime}$. It is observed that for a fixed value of $t_{c}$,
$\langle\Delta{\bf r}^{2}\rangle_{s}$ increases with increase in
$t_{c}^{\prime}$ and for sufficiently large value of $t_{c}^{\prime}$,
$\langle\Delta{\bf r}^{2}\rangle_{s}$ becomes infinitely large, indicating the
escape of the particle without approaching steady state. This observation
suggests that persistence of sufficiently long duration of memory in the
medium provides a kind of elastic bound to the particle that overcomes the
harmonic confinement as a result of which the particle can escape out of the
potential and takes infinitely long time to come back the mean position of the
well or approach the steady state. On the other hand, for a fixed value of
$t_{c}^{\prime}$, $\langle\Delta{\bf r}^{2}\rangle_{s}$ decreases with
increase in $t_{c}$ value and becomes zero for infinitely large value of
$t_{c}$. This fact indicates the trapping of the particle for the presence of
infinitely long duration of activity in the medium. This is the genuine
feature of an active particle confined by a finite potential [49].

Now, substituting the solution $\Xi$ from Eq. (27) and integrating it over all
other variables, we obtain the steady state joint probability distribution
function $P(x,y)$ as

$P(x,y)=\frac{1}{2\pi\sigma^{2}}\exp{\left(-\frac{x^{2}+y^{2}}{2\sigma^{2}}\right)},$
(32)

with the variance $\sigma^{2}$ given by

$\sigma^{2}=\frac{D\left[t_{c}+t_{c}^{\prime}+t_{c}^{2}\Gamma+t_{c}^{\prime
2}(t_{c}+t_{c}^{\prime})\omega_{0}^{2}\right]}{2m^{2}\Gamma\omega_{0}^{2}\left[t_{c}^{\prime}+t_{c}(1+t_{c}\Gamma+t_{c}(t_{c}+t_{c}^{\prime})\omega_{0}^{2})\right]}.$
(33)

In $t_{c}^{\prime}\to 0$ limit, the distribution is still Gaussian with the
variance [Eq. 33] as

$\lim_{t_{c}^{\prime}\to
0}\sigma^{2}=\frac{D(1+t_{c}\Gamma)}{2m^{2}\Gamma\omega_{0}^{2}\left(1+t_{c}\Gamma+t_{c}^{2}\omega_{0}^{2}\right)}.$
(34)

In Figs. 6 (a)-(c), we plot the steady state probability distribution $P(x,y)$
for different values of $t_{c}^{\prime}$ and for a fixed value of $t_{c}=5$.
For all $t_{c}^{\prime}$ values, the distribution is Gaussian centered at the
origin of the potential. As $t_{c}^{\prime}$ increases, $\sigma^{2}$ becomes
an increasing function of $t_{c}^{\prime}$. Hence the variance or width of the
distribution increases and the distribution spreads out [see Figs. 6 (a)-(c)].
At the same time, the peak of the distribution suppresses. This implies that
with increase in memory time scale, the probability of finding the particle at
the mean position of the well decreases and the probability of finding the
particle at larger distances increases. Finally, in $t_{c}^{\prime}\to\infty$
limit, $P(x,y)\to 0$ and the distribution becomes flat [Fig. 6 (c)]. The
vanishing of the position distribution function in $t_{c}^{\prime}\to\infty$
limit represents the escape of the particle out of the potential. This
observation supports the enhancement of particle trajectories [Fig. 3] and the
enhancement of steady state MSD [Fig. 5] with increase of the persistent
duration of memory.

Similarly, in Fig. 6 (d)-(e), we show the plot of $P(x,y)$ for different
values of $t_{c}$ and for a fixed $t_{c}^{\prime}$ value. Initially, in the
$t_{c}\to 0$ limit, $P(x,y$ is Gaussian [Eq. 32] with $\sigma^{2}$ given by

$\lim_{t_{c}\to 0}\sigma^{2}=\frac{D(1+t_{c}^{\prime
2}\omega_{0}^{2})}{2m^{2}\Gamma\omega_{0}^{2}}.$ (35)

As $t_{c}$ is increased, $\sigma^{2}$ becomes decreasing function of $t_{c}$
and hence the variance (or width) of the distribution decreases [see Figs. 6
(d)-(e)]. Simultaneously, the peak of the distribution increases and the
distribution becomes narrow. Finally, in $t_{c}\to\infty$ limit, $\sigma^{2}$
approaches zero value and hence the probability distribution becomes a delta
function since

$\begin{split}\lim_{t_{c}\to\infty}P(x,y)&=\lim_{\sigma\to
0}\frac{1}{2\pi\sigma^{2}}\exp{\left(-\frac{x^{2}+y^{2}}{2\sigma^{2}}\right)}\\\
&=\delta(x)\delta(y).\end{split}$ (36)

This suggests that with increase in $t_{c}$ value, the chances of finding the
particle at the mean position of the well increases, confirming the trapping
of the particle for sufficiently long duration of activity in the medium.
These results are complemented with the suppression of particle trajectories
in Fig. 2 and reduction of steady state MSD in Fig. 5 with increase in $t_{c}$
values.

##  IV SUMMARY AND CONCLUSIONS

In summary, we have explored the self propulsion of a harmonically confined
inertial active Ornstein Uhlenbeck particle in a non-Newtonian environment,
characterized by viscoelastic suspension. We model the dynamics using the
generalized Langevin equation of motion. By solving this non-Markovian model,
we demonstrate a phase diagram distinguishing two separate regimes in
$t_{c}^{\prime}\omega_{0}-t_{c}^{\prime}\Gamma$ parameter space associated
with the oscillatory and non-oscillatory behavior of the solution. From the
exact solution of the dynamics and from the simulation results, both transient
and steady state properties of motion are investigated. The simulated particle
trajectories reveal that the particle exhibits rotational motion in the
oscillatory regime of the parameter space. The rotational trajectories keep on
getting enhanced with increase in persistent duration of memory and gets
suppressed with increase in persistent duration of activity in the medium.
This observation suggests that with increase in memory time scale, the elastic
influence exerted by the environment on the particle dominates over the
harmonic bound of the potential. As a consequence, the particle can come out
of the potential without approaching steady state. This result is further
supported by the enhancement of steady state mean square displacement and
uniform spreading of position distribution function through out the space with
increase in the memory time scale. Similarly, the suppression of particle
trajectories as an increasing function of activity time is also complemented
with the reduction of steady state mean square displacement and narrowing of
the probability distribution function with increase in activity time scale.

##  V Acknowledgement

We thank the 8th statphysics community meeting (ICTS/ISPCM2023/02), during
which some parts of the work were done. MS acknowledges the start-up grant
from UGC, state plan fund from the University of Kerala, and SERB-SURE grant
(SUR/2022/000377) from DST, Govt. of India for financial support. MM
acknowledges SERB international travel grant (ITS/2023/002740) from DST, Govt.
of India for financial support.

## References

  * Bechinger  _et al._ [2016] C. Bechinger, R. Di Leonardo, H. Löwen, C. Reichhardt, G. Volpe, and G. Volpe, Active particles in complex and crowded environments, Rev. Mod. Phys. 88, 045006 (2016).
  * Ramaswamy [2017] S. Ramaswamy, Active matter, J. Stat. Mech. 2017, 054002 (2017).
  * Gompper  _et al._ [2020] G. Gompper, R. G. Winkler, T. Speck, A. Solon, C. Nardini, F. Peruani, H. Löwen, R. Golestanian, U. B. Kaupp, L. Alvarez, T. Kiørboe, E. Lauga, W. C. K. Poon, A. DeSimone, S. Muiños-Landin, A. Fischer, N. A. Söker, F. Cichos, R. Kapral, and e. a. P Gaspard, The 2020 motile active matter roadmap, J. Phys.: Condens. Matter 32, 193001 (2020).
  * Pietzonka [2021] P. Pietzonka, The oddity of active matter, Nat. Phys. 17, 1193 (2021).
  * De Magistris and Marenduzzo [2015] G. De Magistris and D. Marenduzzo, An introduction to the physics of active matter, Physica A 418, 65 (2015).
  * Walther and Müller [2013] A. Walther and A. H. E. Müller, Janus particles: Synthesis, self-assembly, physical properties, and applications, Chemical Reviews 113, 5194 (2013).
  * Howse  _et al._ [2007] J. R. Howse, R. A. Jones, A. J. Ryan, T. Gough, R. Vafabakhsh, and R. Golestanian, Self-motile colloidal particles: from directed propulsion to random walk, Phys. Rev. Lett. 99, 048102 (2007).
  * Wang  _et al._ [2021] G. Wang, T. V. Phan, S. Li, M. Wombacher, J. Qu, Y. Peng, G. Chen, D. I. Goldman, S. A. Levin, R. H. Austin, and L. Liu, Emergent field-driven robot swarm states, Phys. Rev. Lett. 126, 108002 (2021).
  * Lehle and Peinke [2018] B. Lehle and J. Peinke, Analyzing a stochastic process driven by ornstein-uhlenbeck noise, Phys. Rev. E 97, 012113 (2018).
  * Bonilla [2019] L. L. Bonilla, Active ornstein-uhlenbeck particles, Phys. Rev. E 100, 022601 (2019).
  * Martin  _et al._ [2021] D. Martin, J. O’Byrne, M. E. Cates, É. Fodor, C. Nardini, J. Tailleur, and F. van Wijland, Statistical mechanics of active ornstein-uhlenbeck particles, Phys. Rev. E 103, 032607 (2021).
  * Caprini  _et al._ [2019] L. Caprini, U. Marini Bettolo Marconi, A. Puglisi, and A. Vulpiani, Active escape dynamics: The effect of persistence on barrier crossing, The Journal of Chemical Physics 150, 024902 (2019).
  * Caprini and Marini Bettolo Marconi [2021] L. Caprini and U. Marini Bettolo Marconi, Inertial self-propelled particles, J. Chem. Phys. 154, 024902 (2021).
  * Dabelow  _et al._ [2019] L. Dabelow, S. Bo, and R. Eichhorn, Irreversibility in active matter systems: Fluctuation theorem and mutual information, Phys. Rev. X 9, 021009 (2019).
  * Berthier  _et al._ [2019] L. Berthier, E. Flenner, and G. Szamel, Glassy dynamics in dense systems of active particles, The Journal of Chemical Physics 150, 200901 (2019).
  * Wittmann  _et al._ [2018] R. Wittmann, J. M. Brader, A. Sharma, and U. M. B. Marconi, Effective equilibrium states in mixtures of active particles driven by colored noise, Phys. Rev. E 97, 012601 (2018).
  * Fily [2019] Y. Fily, Self-propelled particle in a nonconvex external potential: Persistent limit in one dimension, The Journal of Chemical Physics 150, 174906 (2019).
  * Mandal  _et al._ [2017] D. Mandal, K. Klymko, and M. R. DeWeese, Entropy production and fluctuation theorems for active matter, Phys. Rev. Lett. 119, 258001 (2017).
  * Fodor  _et al._ [2016] É. Fodor, C. Nardini, M. E. Cates, J. Tailleur, P. Visco, and F. van Wijland, How far from equilibrium is active matter?, Phys. Rev. Lett. 117, 038103 (2016).
  * Muhsin  _et al._ [2021] M. Muhsin, M. Sahoo, and A. Saha, Orbital magnetism of an active particle in viscoelastic suspension, Phys. Rev. E 104, 034613 (2021).
  * ten Hagen  _et al._ [2009] B. ten Hagen, S. van Teeffelen, and H. Lowen, Non-gaussian behaviour of a self-propelled particle on a substrate, Condens. Matter Phys. 12, 725 (2009).
  * ten Hagen  _et al._ [2011] B. ten Hagen, S. van Teeffelen, and H. Löwen, Brownian motion of a self-propelled particle, J. Phys.: Condens. Matter 23, 194119 (2011).
  * Cates and Tailleur [2013] M. E. Cates and J. Tailleur, When are active brownian particles and run-and-tumble particles equivalent? consequences for motility-induced phase separation, Euro. Phys. Lett. 101, 20010 (2013).
  * Malakar  _et al._ [2020] K. Malakar, A. Das, A. Kundu, K. V. Kumar, and A. Dhar, Steady state of an active brownian particle in a two-dimensional harmonic trap, Phys. Rev. E 101, 022610 (2020).
  * Buttinoni  _et al._ [2013] I. Buttinoni, J. Bialké, F. Kümmel, H. Löwen, C. Bechinger, and T. Speck, Dynamical clustering and phase separation in suspensions of self-propelled colloidal particles, Phys. Rev. Lett. 110, 238301 (2013).
  * Fily and Marchetti [2012] Y. Fily and M. C. Marchetti, Athermal phase separation of self-propelled particles with no alignment, Phys. Rev. Lett. 108, 235702 (2012).
  * Stenhammar  _et al._ [2014] J. Stenhammar, D. Marenduzzo, R. J. Allen, and M. E. Cates, Phase behaviour of active brownian particles: the role of dimensionality, Soft Matter 10, 1489 (2014).
  * Bialké  _et al._ [2015] J. Bialké, J. T. Siebert, H. Löwen, and T. Speck, Negative interfacial tension in phase-separated active brownian particles, Phys. Rev. Lett. 115, 098301 (2015).
  * Solon  _et al._ [2015] A. P. Solon, J. Stenhammar, R. Wittkowski, M. Kardar, Y. Kafri, M. E. Cates, and J. Tailleur, Pressure and phase equilibria in interacting active brownian spheres, Phys. Rev. Lett. 114, 198301 (2015).
  * Caprini  _et al._ [2021] L. Caprini, C. Maggi, and U. Marini Bettolo Marconi, Collective effects in confined active brownian particles, The Journal of Chemical Physics 154, 244901 (2021).
  * Caprini  _et al._ [2020] L. Caprini, U. M. B. Marconi, C. Maggi, M. Paoluzzi, and A. Puglisi, Hidden velocity ordering in dense suspensions of self-propelled disks, Phys. Rev. Res. 2, 023321 (2020).
  * van Roon  _et al._ [2022] D. M. van Roon, G. Volpe, M. M. Telo da Gama, and N. A. M. Araújo, The role of disorder in the motion of chiral active particles in the presence of obstacles, Soft Matter 18, 6899 (2022).
  * Scholz  _et al._ [2018] C. Scholz, S. Jahanshahi, A. Ldov, and H. Löwen, Inertial delay of self-propelled particles, Nat. Commun. 9, 5156 (2018).
  * Mandal  _et al._ [2019] S. Mandal, B. Liebchen, and H. Löwen, Motility-induced temperature difference in coexisting phases, Phys. Rev. Lett. 123, 228001 (2019).
  * BERG and BROWN [1972] H. C. BERG and D. A. BROWN, Chemotaxis in escherichia coli analysed by three-dimensional tracking, Nature 239, 500 (1972).
  * Martens  _et al._ [2012] K. Martens, L. Angelani, R. Di Leonardo, and L. Bocquet, Probability distributions for the run-and-tumble bacterial dynamics: An analogy to the lorentz model, The European Physical Journal E 35, 84 (2012).
  * Nguyen  _et al._ [2022] G. H. P. Nguyen, R. Wittmann, and H. Löwen, Active Ornstein–Uhlenbeck model for self-propelled particles with inertia, J. Phys.: Condens. Matter 34, 035101 (2022).
  * Noushad  _et al._ [2021] A. Noushad, S. Shajahan, and M. Sahoo, Velocity auto correlation function of a confined Brownian particle, Eur. Phys. J. B 94, 202 (2021).
  * Muhsin and Sahoo [2022] M. Muhsin and M. Sahoo, Inertial active Ornstein-Uhlenbeck particle in the presence of a magnetic field, Phys. Rev. E 106, 014605 (2022).
  * Paraan  _et al._ [2008] F. N. C. Paraan, M. P. Solon, and J. P. Esguerra, Brownian motion of a charged particle driven internally by correlated noise, Phys. Rev. E 77, 022101 (2008).
  * Sevilla  _et al._ [2019] F. J. Sevilla, R. F. Rodríguez, and J. R. Gomez-Solano, Generalized Ornstein-Uhlenbeck model for active motion, Phys. Rev. E 100, 032123 (2019).
  * Gomez-Solano  _et al._ [2022] J. R. Gomez-Solano, R. F. Rodríguez, and E. Salinas-Rodríguez, Nonequilibrium dynamical structure factor of a dilute suspension of active particles in a viscoelastic fluid, Phys. Rev. E 106, 054602 (2022).
  * Narinder  _et al._ [2018] N. Narinder, C. Bechinger, and J. R. Gomez-Solano, Memory-Induced Transition from a Persistent Random Walk to Circular Motion for Achiral Microswimmers, Phys. Rev. Lett. 121, 078003 (2018).
  * Sprenger  _et al._ [2022] A. R. Sprenger, C. Bair, and H. Löwen, Active Brownian motion with memory delay induced by a viscoelastic medium, Phys. Rev. E 105, 044610 (2022).
  * Lozano  _et al._ [2018] C. Lozano, J. R. Gomez-Solano, and C. Bechinger, Run-and-tumble-like motion of active colloids in viscoelastic media, New Journal of Physics 20, 015008 (2018).
  * Narinder  _et al._ [2019] N. Narinder, J. R. Gomez-Solano, and C. Bechinger, Active particles in geometrically confined viscoelastic fluids, New Journal of Physics 21, 093058 (2019).
  * Fox  _et al._ [1988] R. F. Fox, I. R. Gatland, R. Roy, and G. Vemuri, Fast, accurate algorithm for numerical simulation of exponentially correlated colored noise, Phys. Rev. A 38, 5938 (1988).
  * VAN KAMPEN [2007] N. VAN KAMPEN, Chapter viii - the fokker–planck equation, in  _Stochastic Processes in Physics and Chemistry (Third Edition)_, North-Holland Personal Library, edited by N. VAN KAMPEN (Elsevier, Amsterdam, 2007) third edition ed., pp. 193–218.
  * Muhsin and Sahoo [2023] M. Muhsin and M. Sahoo, Inertial active ratchet: Simulation versus theory, Phys. Rev. E 107, 054601 (2023).


# G-SimCLR: Self-Supervised Contrastive Learning with Guided Projection via
Pseudo Labelling

Souradip Chakraborty§ Walmart Labs
<EMAIL_ADDRESS>
Aritra Roy Gosthipaty§ Netaji Subhash Engineering College
<EMAIL_ADDRESS>Sayak Paul§ PyImageSearch
<EMAIL_ADDRESS>
###### Abstract

In the realms of computer vision, it is evident that deep neural networks
perform better in a supervised setting with a large amount of labeled data.
The representations learned with supervision are not only of high quality but
also helps the model in enhancing its accuracy. However, the collection and
annotation of a large dataset are costly and time-consuming. To avoid the
same, there has been a lot of research going on in the field of unsupervised
visual representation learning especially in a self-supervised setting.
Amongst the recent advancements in self-supervised methods for visual
recognition, in SimCLR Chen et al. shows that good quality representations can
indeed be learned without explicit supervision. In SimCLR, the authors
maximize the similarity of augmentations of the same image and minimize the
similarity of augmentations of different images. A linear classifier trained
with the representations learned using this approach yields 76.5% top-1
accuracy on the ImageNet ILSVRC-2012 dataset. In this work, we propose that,
with the normalized temperature-scaled cross-entropy (NT-Xent) loss function
(as used in SimCLR), it is beneficial to not have images of the same category
in the same batch. In an unsupervised setting, the information of images
pertaining to the same category is missing. We use the latent space
representation of a denoising autoencoder trained on the unlabeled dataset and
cluster them with $k$-means to obtain pseudo labels. With this apriori
information we batch images, where no two images from the same category are to
be found. We report comparable performance enhancements on the CIFAR10 dataset
and a subset of the ImageNet dataset
111https://github.com/thunderInfy/imagenet-5-categories. We refer to our
method as G-SimCLR 222Code available at https://github.com/ariG23498/G-SimCLR
333This paper is accepted at ICDM 2020 DLKT workshop..

###### Index Terms:

Self-supervised Representation Learning, Clustering, Contrastive Learning,
Visual Recognition

§§footnotetext: Equal contribution

##  I Introduction

The huge success of deep learning algorithms in computer vision is heavily
dependent on the quality and the information content of the latent visual
representations learned from the input images. There has been almost a decade
of very successful research in the field of supervised representation learning
from images using different variants of convolutional neural networks [1, 2,
3, 4, 5]. However, the complex parametric structure of deep convolutional
neural networks makes it heavily dependent on the amount of labeled training
data. Though, humongous amount of labeled data helps in increasing both the
accuracy of these models and the information content of the visual
representation learned in these models, the collection and labeling of a large
dataset are extremely expensive and time-consuming. However, the amount of
unlabeled data (e.g. images on the Internet) is practically unlimited in
nature except for some domains like healthcare. Hence to remedy the expensive
efforts of creating a huge amount of labeled datasets, there has been an
increasing amount of research dedicated to the field of unsupervised (i.e.
where no explicit label information is used) visual representation learning.
Methods that have come out as a byproduct of this research effort, have shown
tremendous performance with respect to their fully supervised counterparts.
Features learned by these methods have been very useful for several downstream
visual recognition tasks such as image classification, object detection,
instance segmentation, and so on even with very less labeled data [6, 7, 8, 9,
10].

When dealing with unsupervised regimes, methods that are based on self-
supervised learning have shown great promise. In self-supervised learning,
explicit labels are not used; instead, we utilize the unlabeled dataset to
formulate an auxiliary supervised learning problem e.g. predicting rotation
angle from images, placing image patches in the right place, placing video
frames in the right order, and so on [11, 12, 13].

Another approach that we see in SimCLR [9] is contrasting different views of
images and training models to maximize the agreement between the views of the
same images while minimizing the agreement between the views of other images.
A linear classifier trained with the representations learned using this
approach yields a 76.5% top-1 accuracy on the ImageNet ILSVRC-2012 dataset
[14].

However, in [15] Khosla et al. show that SimCLR being dependent on the
sampling probability might learn representations that can distant images from
the same label space if they appear in the same training batch. This in turn
hurts the representations learned for subsequent downstream tasks. This can
turn out to be a crucial problem if the downstream task is highly imbalanced
in its label space. For example, in a dataset with 5 classes if one class is a
majority class having the maximum number of training samples, then there is a
very high probability that samples from that majority class will be in the
same batch while training SimCLR. If images belonging to the same class are
sent to the same batch for training in SimCLR, the model might learn
representations that can distant the images of the same class in the latent
space which can hurt the quality of the learned representations.

Motivated by these findings and work present in [16] by Caron et al. we
propose a simplistic approach for generating prior information of the label
space without using any explicit labels and then trying to enhance the quality
of the representations learned via SimCLR through this apriori information. We
use a denoising autoencoder (fully convolutional) [17] in a self-supervised
setting where we minimize the reconstruction loss between the input and noisy
images and extract the latent space representations from the trained
autoencoder. We leverage this latent feature map from the autoencoder to
cluster the input space of images and assign each input image to the
corresponding cluster which we refer to as pseudo labels. We incorporate this
prior pseudo label information while determining the training batches for
SimCLR which minimizes the risk of getting images from the same label space in
the same batch. We report comparable performance gains achieved using G-SimCLR
i.e. the proposed method on two datasets - CIFAR10 [18] and a subset of
ImageNet [19] (referred to as ImageNet Subset throughout the paper). To
encourage research along similar lines we also open-source the code to
reproduce our experiments.

##  II Methodology

In SimCLR, Chen et al. have shown that a strong suite of stochastic data
augmentation operations, suitable non-linear transformations and contrastive
learning with large batch size, and longer training time enhance the quality
of the learned representations significantly. In this work, we present an
additional step on top of the SimCLR framework to ensure two semantically
similar images do not get treated differently. The schematics of this
additional step is presented below -

  * •

First, we train a denoising autoencoder (fully convolutional) on the given
dataset.

  * •

In the second step, we take the encoder-learned robust representations (also
referred to as latent space representations) and use $k$-means to get initial
cluster assignments of those representations. We refer to these cluster
assignments as pseudo labels for the given dataset.

  * •

Finally, we use these pseudo labels to prepare the batches and use them for
SimCLR training.

We present this methodology in Figure 1. In the first step of our methodology,
a denoising autoencoder is used which enhances the robustness of the latent
feature representation when compared to a vanilla autoencoder. This is done by
distorting the input image(x) to get a semi-distorted version of the same by
means of Gaussian noise. The corrupted input($\hat{x}$) is then mapped, as in
case of a vanilla autoencoder, to a latent representation $z=f(\hat{x})$,
where $z\in R^{d}$. The latent encoder representation $z$ now helps in
reconstructing the distortion-free image input in the decoding phase and the
average reconstruction loss is minimized to learn the parameters of the
denoising autoencoder.

Once the denoising autoencoder is trained, the encoder representations are
extracted from all the input images. This representation of the input data
space is denoted by $Y\in R^{n\times d}$, where $n$ is the total number of
input images and $d$ is the dimensionality of the encoder representation. We
use $k$-means on this representation space of the images ($Y$) to get dense
clusters (assume the number of clusters to be $k$). The clusters formed by
minimizing the within-group variance and maximizing the between-group
variance, ensures that the images which are similar to each other fall in the
same cluster and which are different falls in different clusters. The cluster
assignment is given by $d_{ij}=1$, when $i_{th}$ image falls in the $j_{th}$
cluster else $d_{ij}=0$. These assignments are considered as pseudo labels in
our work, used for selecting the images for each batch in training SimCLR.

In the final step of our methodology, instead of randomly sampling the mini-
batch of $p$ samples, we use the pseudo labels (obtained by $k$-means
clustering on the autoencoder representations) to sample the images in each of
the batches. In order to minimize the risk of having similar images in the
same batch while SimCLR training, which might adversely impact the
representations, we use a pseudo label based stratified sampling methodology
in preparing the training batches. In our current implementation, we have kept
the batch size and cluster number to be the same, $k$. So, while preparing the
training batches, we perform random sampling (without replacement), stratified
by the pseudo labels. $k$-means clustering does not ensure equal-sized
clusters. This means that some batches will have more from one cluster and
less from others. This ensures that our method does not provide a hard
constraint on the images of the batches to belong to different and discrete
classes only.

We then follow the same training objective as SimCLR but on the pairs of
augmented examples derived from the batches formed using the aforementioned
methodology. The loss function (NT-XEnt) for training SimCLR (as well as
G-SimCLR) is expressed as

$\ell_{i,j}=-\log\frac{\exp\left(\operatorname{sim}(\boldsymbol{z}_{i},\boldsymbol{z}_{j})/\tau\right)}{\sum^{2N}_{k=1}\mathbbm{1}_{[k\neq
i]}\exp\left(\operatorname{sim}(\boldsymbol{z}_{i},\boldsymbol{z}_{k})/\tau\right)}$
(1)

where $sim$ denotes cosine similarity, $z$ denotes the outputs from the non-
linear projection head as used in the original SimCLR work, and $\tau$ denotes
the temperature hyperparameter.

Overall, the primary additional step as discussed above does introduce an
additional bottleneck in terms of computational complexity and execution time.
We do not put too much focus on generating the pseudo labels, these are used
only to guide the batch formation systematically. In the original SimCLR work
large batch sizes (typically 4096) are used to have enough negative samples.
However, maintaining such large batch sizes is computationally expensive and
G-SimCLR provides a simple way to remedy that.

Figure 1: The proposed method involves training a denoising autoencoder to
capture a latent space representation for the given dataset. We then cluster
this latent space to generate pseudo labels for the given dataset. Finally, we
use these pseudo labels to prepare the data batches for SimCLR training on the
given dataset.

##  III Experimental setup

Denoising autoencoder. A Gaussian noise centered at 0 with a standard
deviation of 0.01 is provided as the input layer of the denoising autoencoder.
This layer adds Gaussian noise to the input images. The autoencoder learns
from the mean squared error loss (MSE loss) between the noisy and the original
images. Adam optimizer [20] is used to update the weights of the network.
Using early stopping (with a patience of 5) the autoencoder trains for 100
epochs. After the autoencoder is trained, the latent space representation of
each training image is extracted. The representations are taken from the
bottleneck of the autoencoder.

Clustering with $k$-means. By applying $k$-means clustering algorithm on the
latent space representations of images, we retrieve the cluster labels for
each image. The choice for the number of clusters is 64. Empirically, we found
out that it is helpful to have a cluster number that is large enough to
capture the implicit marginal distribution of the given dataset. One could
also use domain knowledge to determine this value [21]. More rigorous
experimentation to study the effect of the number of $k$ on the contrastive
learning task and further downstream tasks has been left as a scope for future
work.

SimCLR and G-SimCLR. Our approach of SimCLR training is not exactly similar to
what is presented in the original paper. We refer to this variant of SimCLR as
SimCLR with minor modifications. ResNet [4] is used as the feature backbone
network to extract meaningful representations from the input images (and their
augmentation variants). Similar to SimCLR we map these representations to a
lower-dimensional space using a shallow fully-connected network with ReLU non-
linearity. We use a subset of augmentation operations used in SimCLR.
Specifically, we use random horizontal flips and color distortions. The
optimizer used is stochastic gradient descent with a cosine learning rate
decay [22] instead of LARS [23]. A temperature ($\tau$) value of 0.1 is used.

Linear evaluation and label efficient fine-tuning. For linear evaluation of
the representations learned using our variant of SimCLR and G-SimCLR we only
use a linear layer. We further fine-tune the representations with 10% labeled
examples from the datasets we use. We report the final performances from these
experiments in Tables II and III respectively. For these experiments, we use
early stopping and train with the default values of the Adam optimizer in
Keras 444https://keras.io/api/optimizers/adam/.

Our variant of SimCLR and G-SimCLR differ only in the way the training batches
are made. In our proposed approach, the training images are ordered in a batch
of 64 samples. The organization of images in a batch is such that no two
images belonging to the same cluster are to be found in the same batch. The
choice of 64 as the batch size is not cherry-picked. We use 64 clusters so
that in a batch of 64 images, we can group images belonging to distinct
clusters in a batch.

Below we present dataset-specific configuration details.

###  III-A CIFAR10

The latent space of the autoencoder is (4,4,128) dimensional. Our variant of
SimCLR and G-SimCLR are trained for 15 epochs. For the feature backbone
network, we use a variant of ResNet that is specifically suited for the
CIFAR10 dataset 555We use this implementation -
https://github.com/GoogleCloudPlatform/keras-idiomatic-
programmer/blob/master/zoo/resnet/resnet_cifar10.py.

###  III-B ImageNet Subset

We experiment with two autoencoders for this dataset. With the first kind, we
use 3 hidden layers for encoder and decoder respectively. The latent space
here is (28,28,128) dimensional. The second kind has 5 hidden layers in the
encoder and the decoder. The latent space here is (7,7,1024) dimensional. We
find that the deeper autoencoder encapsulates more robust and discriminative
visual representations in its latent space. We use ResNet50 666We use
tf.keras.applications.ResNet50 as seen here -
https://keras.io/api/applications/resnet/. as the feature backbone network and
we train it for 200 epochs.

##  IV Experimental results

We observe the performance of the fully convolutional denoising autoencoder on
a given dataset in terms of its reconstruction loss. Since our approach is
completely unsupervised in nature, the train and validation reconstruction
losses give us a directional understanding of the quality of the
representations learned with the autoencoder network as shown in 2. The
convergence of the train and validation reconstruction losses indicate that
the representations learned are meaningful and informative.

Figure 2: Train and validation reconstruction losses of the denoising
autoencoder networks trained on the CIFAR10 and ImageNet Subset datasets
respetively.

In an unsupervised setting, we do not know how many categories the images
belong to. In our approach, we cluster the latent space representations of
images into 64 clusters. The image representations close to each other fall
into the same cluster, while the ones that are far away, belong to different
clusters. From the t-SNE [24] projections (as shown in Figure 3) of the latent
space representations (as learned by the denoising autoencoder), it is evident
that the representations of similar images cluster together. We do not make
use of the typical cluster evaluation metrics like Normalized Mutual
Information Score [25] and Rand Index [26] since they involve comparisons with
the ground-truth labels, thereby eliminating the need for supervision as much
as possible. We do not treat the clustering problem very rigorously; we make
use of it just enough to generate the pseudo labels needed to systematically
construct the mini-batches for SimCLR.

Figure 3: t-SNE projections of the latent space representations learned by the
denoising autoencoder on the CIFAR10 and ImageNet Subset datasets.

In Table I we show how we represent the pseudo labels for G-SimCLR.

TABLE I: Image indices and their cluster labels as generated by $k$-means. We treat these cluster labels as pseudo labels in G-SimCLR. Image index | Cluster label
---|---
36 | 0
61 | 15
99 | 47
11 | 48
47 | 14

Figure 4 shows the loss (NT-Xent) curves as obtained from the G-SimCLR
training with the CIFAR10 and ImageNet Subset datasets. We can see that in
both cases the loss decrease is steady suggesting good progress in G-SimCLR
training.

Figure 4: Loss (NT-Xent) curves as obtained from the G-SimCLR training with
the CIFAR10 and ImageNet Subset datasets respectively.

To compare G-SimCLR with our variant of SimCLR, we run two downstream
supervised classification tasks on the visual representations that each of
them learns. We first run linear evaluation on the classification task and
later fine-tune the models with 10 % data and then evaluate again.

For linear evaluation we report the validation accuracy in Table II for
different levels of the feature backbone network and the projection head
denoted as P1, P2 and P3 where

  * •

P1 denotes the feature backbone network + the entire non-linear projection
head - its final layer

  * •

P2 denotes the feature backbone network + the entire non-linear projection
head - its final two layers

  * •

P3 denotes the feature backbone network only

Next in Table III we report the validation accuracy where the representations
learned by our variant of SimCLR and G-SimCLR are fine-tuned with only 10%
labeled training examples. In this case, we do not use any layer from the non-
linear projection head. The results reported suggest that the intuition behind
G-SimCLR is valid.

The results reported suggest that the guided apriori provided to batches in
G-SimCLR is indeed more performant than our SimCLR variant.

TABLE II: Performance of the linear classifiers trained on top of the
representations (kept frozen during training linear classifiers) learned by
G-SimCLR. Linear evaluation
---
|  | CIFAR10 | ImageNet Subset
Fully supervised |  | 73.62 | 67.6
SimCLR with minor modifications | P1 | 37.69 | 52.8
P2 | 39.4 | 48.4
P3 | 39.92 | 52.4
G-SimCLR (ours) | P1 | 38.15 | 56.4
P2 | 41.01 | 56.8
P3 | 40.5 | 60.0
TABLE III: Performance of weakly supervised classifiers trained on top of the
representations learned by G-SimCLR. We fine-tuned the representations with
only 10% labeled data. Fine-tuning (10% labeled data)
---
| CIFAR 10 | ImageNet Subset
SimCLR with minor modifications | 42.21 | 49.2
G-SimCLR (ours) | 43.1 | 56.0

##  V Conclusion

In this work, we build on top of the original SimCLR work and present a simple
yet effective way to form robust mini-batches for SimCLR training. We study
the main components used in G-SimCLR (proposed method) i.e. denoising
autoencoder and the intermediary clustering step and we report comparable
performance enhancements on two datasets.

We are hopeful G-SimCLR will be able to remedy the problems of SimCLR e.g.
requirement of maintaining large batch sizes, false positives in the learned
representation space, and so on. Benchmarking G-SimCLR on more standardized
datasets (such as ImageNet, Flowers, Caltech-101, etc.), more detailed study
on the autoencoder-learned latent space representations, using more efficient
algorithms for clustering remain as future scopes of improvements for us.

## References

  * [1] A. Krizhevsky, I. Sutskever, and G. E. Hinton, “Imagenet classification with deep convolutional neural networks,” in Advances in Neural Information Processing Systems 25 (F. Pereira, C. J. C. Burges, L. Bottou, and K. Q. Weinberger, eds.), pp. 1097–1105, Curran Associates, Inc., 2012.
  * [2] K. Simonyan and A. Zisserman, “Very deep convolutional networks for large-scale image recognition,” 2014.
  * [3] C. Szegedy, V. Vanhoucke, S. Ioffe, J. Shlens, and Z. Wojna, “Rethinking the inception architecture for computer vision,” 2015.
  * [4] K. He, X. Zhang, S. Ren, and J. Sun, “Deep residual learning for image recognition,” 2015.
  * [5] M. Tan and Q. V. Le, “Efficientnet: Rethinking model scaling for convolutional neural networks,” 2019.
  * [6] T. H. Trinh, M.-T. Luong, and Q. V. Le, “Selfie: Self-supervised pretraining for image embedding,” 2019.
  * [7] I. Misra and L. van der Maaten, “Self-supervised learning of pretext-invariant representations,” 2019.
  * [8] K. He, H. Fan, Y. Wu, S. Xie, and R. Girshick, “Momentum contrast for unsupervised visual representation learning,” 2019.
  * [9] T. Chen, S. Kornblith, M. Norouzi, and G. Hinton, “A simple framework for contrastive learning of visual representations,” 2020.
  * [10] J.-B. Grill, F. Strub, F. Altché, C. Tallec, P. H. Richemond, E. Buchatskaya, C. Doersch, B. A. Pires, Z. D. Guo, M. G. Azar, B. Piot, K. Kavukcuoglu, R. Munos, and M. Valko, “Bootstrap your own latent: A new approach to self-supervised learning,” 2020.
  * [11] S. Gidaris, P. Singh, and N. Komodakis, “Unsupervised representation learning by predicting image rotations,” 2018.
  * [12] M. Noroozi and P. Favaro, “Unsupervised learning of visual representations by solving jigsaw puzzles,” 2016.
  * [13] H.-Y. Lee, J.-B. Huang, M. Singh, and M.-H. Yang, “Unsupervised representation learning by sorting sequences,” 2017.
  * [14] O. Russakovsky, J. Deng, H. Su, J. Krause, S. Satheesh, S. Ma, Z. Huang, A. Karpathy, A. Khosla, M. Bernstein, A. C. Berg, and L. Fei-Fei, “ImageNet Large Scale Visual Recognition Challenge,” International Journal of Computer Vision (IJCV), vol. 115, no. 3, pp. 211–252, 2015.
  * [15] P. Khosla, P. Teterwak, C. Wang, A. Sarna, Y. Tian, P. Isola, A. Maschinot, C. Liu, and D. Krishnan, “Supervised contrastive learning,” 2020.
  * [16] M. Caron, P. Bojanowski, A. Joulin, and M. Douze, “Deep clustering for unsupervised learning of visual features,” 2018.
  * [17] P. Vincent, H. Larochelle, Y. Bengio, and P.-A. Manzagol, “Extracting and composing robust features with denoising autoencoders,” in International Conference on Machine Learning proceedings, 2008.
  * [18] A. Krizhevsky, “Learning multiple layers of features from tiny images.” https://www.cs.toronto.edu/~kriz/learning-features-2009-TR.pdf, 2009.
  * [19] A. Rastogi, “Imagenet 5 categories.” https://github.com/thunderInfy/imagenet-5-categories, 2020.
  * [20] D. P. Kingma and J. Ba, “Adam: A method for stochastic optimization,” 2014.
  * [21] Y. M. Asano, C. Rupprecht, and A. Vedaldi, “Self-labelling via simultaneous clustering and representation learning,” 2019.
  * [22] I. Loshchilov and F. Hutter, “Sgdr: Stochastic gradient descent with warm restarts,” 2016.
  * [23] Y. You, I. Gitman, and B. Ginsburg, “Large batch training of convolutional networks,” 2017.
  * [24] L. van der Maaten and G. Hinton, “Visualizing data using t-SNE,” Journal of Machine Learning Research, vol. 9, pp. 2579–2605, 2008.
  * [25] A. Strehl and J. Ghosh, “Cluster ensembles — a knowledge reuse framework for combining multiple partitions,” JMLR, 2003.
  * [26] W. Rand, “Objective criteria for the evaluation of clustering methods. journal of the american statistical association.,” JASA, 1971.