# Effect of fluctuations on the Geodesic rule for topological defect formation

Sanatan Digal<EMAIL_ADDRESS>The Institute of Mathematical Sciences,
Chennai 600113, India  Homi Bhabha National Institute, Training School
Complex, Anushakti Nagar, Mumbai 400085, India  Vinod Mamale
<EMAIL_ADDRESS>The Institute of Mathematical Sciences, Chennai 600113,
India  Homi Bhabha National Institute, Training School Complex, Anushakti
Nagar, Mumbai 400085, India

###### Abstract

At finite temperature, the field along a linear stretch of correlation length
size is supposed to trace shortest path in the field space given the two end
points, known as the Geodesic rule. In this study we compute the probability
that the field variations over distances of correlation length follow the
Geodesic rule in theories with $O(2)$ global symmetry. We consider a model of
ferromagnetic $O(2)$ spins and a complex $\phi^{4}$ theory. The computations
are carried out on an ensemble of equilibrium configurations at finite
temperatures generated using Monte Carlo simulations. The numerical results
suggest that, for temperatures relevant for the studies of topological defect
formation during 2nd order phase transition, there is a significant deviation
from the Geodesic rule. We also study the equilibrium density and distribution
of vortices in O(2) spins in two dimensions and compare with the Kibble-Zurek
Mechanism of defect formation.

##  I Introduction

Topological defects arise in a wide variety of systems ranging from table-top
experiments in condensed matter systemsChuang:1991zz ; Hendry:1994 ;
Ruutu:1996 ; Bauerle:1996 ; Dodd:1998aan ; Carmi:2000zz ; Digal:1998ak  to
theories of the early Universe, extremely dense stars etc.Hindmarsh:1994re ;
Magueijo:2000se ; Srivastava:2017itj . Topological defects mostly form during
phase transitions accompanied by spontaneous breaking (SSB) of discrete or
continuous symmetries. These defects can also be excited in the SSB phase by
external stimulations Digal:1999ub . The type of topological defects that are
possible depends on the corresponding homotopy group and the spatial
dimensions, irrespective of the energetics involved Zurek:1996sj ;
Toulouse:1976 ; Mermin:1979zz . The processes of formation and evolution of
these defects have been extensively studied in the literature. There are many
theoretical as well as experimental studies available on the formation of
topological defects in condensed matter systems Srivastava:2001 ; Chaikin2000
. In the case of high energy physic systems the studies are mostly
theoretical. However since topological considerations play a dominant role in
their formation and evolution, ideas proposed for high energy physics systems
have been tested in condensed matter systems.

The theory of formation of topological defects was first proposed in the
context of the early Universe by T. W. B. Kibble, known as the Kibble-
Mechanism(KM) Kibble:1976sj . The Kibble-Mechanism is based on two postulates.
In the immediate aftermath of a phase transition the order parameter (OP) in
physical space takes values from the order parameter space (OPS). According to
the first postulate of KM, immediately after a SSB phase transition the
physical space splits into domains, inside of which the OP is roughly uniform.
The second postulate, which is also known as the “Geodesic Rule”, states that
in between adjacent domains OP field interpolates along the shortest path in
the OPS Kibble:1995aa . Given an OPS these two postulates can be applied to
make definite predictions such as the minimum number of domains required to
form a defect, defect densities and defect correlations etc. Bowick:1992rz ;
Ray:2001ug . For example, in two dimensions a minimum of three domains are
required to form a vortex or anti-vortex, which is usually located near the
junction of these. The probability of formation is estimated to be $0.25$.
According to the Kibble Mechanism the formation of defects for a given OPS,
spatial dimensions, and the nature of phase transition, does not depend on the
energy scales involved, i.e whether in a condensed matter system, in high
energy systems or in the early Universe.

The Kibble mechanism has been very successful in describing the formation of
topological defects during a first order transition. In a first order phase
transition, bubbles of SSB phase nucleate in the background of symmetric phase
Srivastava:1991nv . The bubbles with uniform OP always dominate during the
nucleation as these correspond to least action. These bubbles can be
identified with the domains proposed in the Kibble Mechanism. These bubbles
grow, collide/coalesce with one another converting the symmetric phase into
the SSB phase. It has been observed in numerical experiments that when two
bubbles collide, the OP variation from one bubble to the other interpolates
along the shortest path in the OPS validating the geodesic rule of the Kibble
mechanism Srivastava:1991nv ; Chakravarty:1992zp . Though there can be
exceptions to this when kinematics allow energetic bubble collisions or there
is a small explicit symmetry breaking present apart from the SSB
Copeland:1996jz ; Digal:1997ip ; Digal:1995vd .

In a second order phase transition formation of topological defects is
described by the Kibble-Zurek mechanism Zurek:1996sj . In this case it becomes
necessary to incorporate the phenomenon of critical slowing down in the Kibble
Mechanism which was first pointed out by Zurek. Zurek argued that in a second
order phase transition as the system temperature approaches the critical
temperature$(T_{c})$ from above, the dynamics of the field almost freezes as
soon as it enters the regime of critical slowing down
$\left[T_{c}+\epsilon_{1},T_{c}-\epsilon_{2}\right]$. When the system exits
this regime, the field configuration still corresponds to temperature
$T_{c}+\epsilon_{1}$, even though the temperature is equal to or below
$T_{c}-\epsilon_{2}$. The field configuration at $T=T_{c}+\epsilon_{1}$ is
laid out using the two postulates of Kibble mechanism. If the temperature of
the system were to remain fixed at $T_{c}+\epsilon_{1}$ the domains will keep
fluctuating and there will be continuous creation and anihilation of unstable
topological defects. Since the system temperature $T\leq T_{c}-\epsilon_{2}$,
i.e the effective potential now has a non trivial OPS, the field inside the
domains will roll down to the corresponding nearest point on the OPS. In this
process the topological defects will become stable.

Though the two postulates of the Kibble Mechanism and the Kibble-Zurek
Mechanism, have been very successful in describing the formation and dynamics
of topological defects their validity, in particular that of the Geodesic rule
has never been tested in the case of $2nd$ order phase transition. It is
argued that the Geodesic rule follows from free energy considerations, i.e
when the field variation over the distances traces shortest path on the OPS
the energy and free energy are minimised. However, the same can be achieved by
maximising the entropy. A definitive answer to whether there are deviations to
Geodesic rule can come only from the simulations of the partition function. In
this work we compute the deviations from the Geodesic rule in the $O(2)-$spin
model in two and three spatial dimensions. In this model the spins take values
from an unit circle. We also consider $\phi^{4}$ theory with $U(1)\equiv O(2)$
symmetry in three spatial dimensions to see the effect of radial/magnitude
fluctuations on the Geodesic rule. In the relevant temperature range of
interest the $\phi$ field takes values from a circular band in the complex
$\phi-$plane.

The computations involve Monte Carlo simulation of the lattice partition
function of both the models. At a given temperature the correlation length
$\xi$ is computed. The conventional definition of $\xi$ is not suitable for
the Kibble Mechanism as the correlation of the field at $\xi$ separations does
not vanish. Therefore we define $\xi_{KM}$ such that field correlation
approximately vanishes over $\xi_{KM}$ separations. It is found that
$\xi_{KM}$ is significantly larger than $\xi$. Once $\xi$ and $\xi_{KM}$ are
computed the field configuration over a stretch of $\xi_{KM}$ is mapped onto
the field space. Finding length of the image trajectory in the field space
requires assumption of the Geodesic rule for variations of the field between
nearest neighbour (NN) on the lattice. The simulations results for both the
models show significant deviations from the Geodesic rule. The deviation over
$\xi$ distances is found to be above $20\%$. Over $\xi_{KM}$ the trajectories
no longer prefer the Geodesic path in the field space. Even the maximal
violation of the Geodesic rule largely underestimates the defect densities.
This is possibly due to dominance of defects arising from thermal
fluctuations. These results suggest that even over $\xi_{KM}$ fluctuations
dominate the field dynamics.

The paper is organised as follows. In section-II and section-III, numerical
calculations to determine the validity of geodesic rule are presented for the
$O(2)-$spin model and $\phi^{4}$ theory respectively. The conclusions are
presented in section-IV.

##  II Geodesic rule in $O(2)-$spin model

The $O(2)-$spin model, also known as the $XY-$model, is a toy spin model on
the lattice Kosterlitz:1973xp  . Each of the spins, $\vec{s}_{i}$ at the site
$i$, is a two component unit vector. The spin $\vec{s}_{i}$ can also be
represented by an angle $\theta_{i}$ such that
$\vec{s}_{i}=\left(cos\theta_{i},sin\theta_{i}\right)$. The Hamiltonian,

$H=-J\sum_{\left<ij\right>}\vec{s}_{i}.\vec{s}_{j},$ (1)

includes ferromagnetic coupling ($J>0$) between the $i-$th and $j-$th spins
which are NN’s. The $O(2)-$symmetry corresponds to global rotations of the
spins which preserve the Hamiltonian. For simplicity, external field breaking
the $O(2)$ symmetry explicitly is not considered in the present study. The
expression for the partition function is given by,

${\cal{Z}}=\int\prod_{i}^{N}d\vec{s}_{i}Exp[\beta\sum_{\left<ij\right>}\vec{s}_{i}.\vec{s}_{j}].$
(2)

$\beta=J/\kappa T$, where $\kappa$ and $T$ are the Boltzmann constant and
temperature respectively. For convenience $T$ is considered in units of
$J/\kappa$. It is well known that in three or higher dimensions the system of
lattice spins undergoes a $2$nd order ferromagnetic-paramagnetic transition at
critical temperature $T_{c}$. The order parameter for this transition is given
by,

$\vec{M}=\left<{1\over N}\sum_{i}\vec{s}_{i}\right>.$ (3)

The magnitude of $\vec{M}$ acquires non-zero value for temperatures, $T\leq
T_{c}$, which leads to SSB of $O(2)$ symmetry. The symmetry broken phase
allows for existence of stable topological defects corresponding to the first
homotopy group $\pi_{1}(S^{1})\equiv Z$. These defects can exist even above
$T_{c}$ but are unstable. In the case of two physical dimensions there is no
magnetisation transition due to dominance of the Goldstone modes.
Interestingly the system undergoes the well known Kosterliz-Thouless(KT)
transitionKosterlitz:1973xp . The spin system consists of a network of
vortices at all temperatures. The two phases are characterised by how the
vortices interact. The high temperature phase is described by a screened
Coulomb gas and a strongly interacting system below the KT transition. Though
the two dimensional system is not directly relevant for the formation of
topological defects via the Kibble-Zurek mechanism, the distribution of
vortices at high temperature serves as an instructive example. The results of
the deviation from the Geodesic rule in two and three dimensions are very
similar. Also there is no reason to expect that the Geodesic rule will not
hold in the case of two spatial dimensions.

Simulation of the partition function, Eq.2, involves mainly generating
statistically significant spin configurations using Monte Carlo methods at
various temperatures above $T>T_{c}$. To generate the configurations the
standarad metropolis algorithmHasting:1970  is used. In this algorithm the
spin $\vec{s}_{i}$ at site $i$ is rotated by an angle $\theta$ randomly chosen
between $0-2\pi$ to obtain $\vec{s}^{\prime}_{i}$. The move,
$\vec{s}_{i}\rightarrow\vec{s}_{i}^{\prime}$ is accepted if the change in the
energy $\Delta E=E^{\prime}-E<0$. If $\Delta E>0$ then the move is accepted
with probability $Exp(-\beta\Delta E)$. While the spin $\vec{s}_{i}$ is being
updated, the rest of the spins are held fixed. A sweep constitutes the
sequential updating of all the spins. Since a new spin configurations is
obtained from an old one, there is non-zero autocorrelation between them. To
reduce this auto-correlation to a reasonable level a configuration is accepted
for computations after about $~{}10$ sweeps. For simplicity, a square lattices
($N=N_{s}\times N_{s}$) and a cubic lattice ($N=N_{s}\times N_{s}\times
N_{s}$) are considered in $2$ and $3-$ dimensions respectively. In Fig.2 and
Fig.2 show a small part of typical 2-dimensional spin configuration generated
at temperature $T=1.43$ and $T=1.25$ respectively. In these figures the purple
dots represent vortices and green dots represent anti-vortices. The spin
configuration shown in Fig.2 has more fluctuations compared to configuration
in Fig.2 as the temperature is higher for the former. The larger fluctuations
lead to higher defect density, so it grows with temperature.

Figure 1: A sampled spin configuration
at $T=1.6T_{c}$.

Figure 2: A sampled spin configuration
at $T=1.4T_{c}$.

###  II.1 Calculation of $\xi$ and $\xi_{KM}$

Conventionally the correlation length $\xi$ is obtained from the correlation
function,

$c(r)=\left\langle\vec{s}_{i}\cdot\vec{s}_{j}\right\rangle-\left\langle\vec{s}_{i}\right\rangle\cdot\langle\vec{s}_{j}\rangle,~{}~{}r=|{\rm\bf
r}_{i}-{\rm\bf r}_{j}|,$ (4)

of spins $\vec{s}_{i}$ and $\vec{s}_{j}$ located at ${\rm\bf r}_{i}$ and
${\rm\bf r}_{j}$ respectively. For large $r$, $c(r)$ is expected to decay as
$\sim Exp[-r/\xi(T)]$. Clearly at $r=\xi$, $c(r)\sim 1/e$, consequently
$\vec{s}_{i}$ can not take values independent of $\vec{s}_{j}$ and vice-versa.
For $T=1.18$, $\xi\simeq 4.4$ in lattice units. At $\xi$ separations, the
probability that the relative angle
$\Delta\theta=\rm{min}\left[|\theta_{i}-\theta_{j}|,2\pi-|\theta_{i}-\theta_{j}|\right]$
between $\vec{s}_{i}$ and $\vec{s}_{j}$ is below $\pi/2$ is $\sim 65\%$, i.e
$\langle\Delta\theta\rangle\neq\pi/2$ and
$\langle\theta_{i}\theta_{j}\rangle\neq\pi^{2}$. According to the Kibble
Mechanism $\Delta\theta$, between two adjacent domains, can take any value
between $0-\pi$ with uniform probability. Consequently,
$\langle\Delta\theta\rangle=\pi/2$ or
$\langle\theta_{i}\theta_{j}\rangle=\pi^{2}$. Therefore $\xi$ does not
accurately describe the size of the domains prescribed in the Kibble
Mechanism. We define $\xi_{KM}$ as the least physical separation at which the
two spins $\vec{s}_{i}$ and $\vec{s}_{j}$ are not correlated, i.e
$c(\xi_{KM})\simeq 0$. At $T=1.18$, it so happens that $\xi_{KM}$ is larger
than $\xi$ by a factor of $\sim 5$.

###  II.2 The Geodesic rule over $\xi$ and $\xi_{KM}$ separations

The calculations to check the validity or the deviations of the Geodesic rule
necessitates that given a configuration of the field along a linear stretch in
physical space is mapped to a trajectory in the field space. Since the field
is defined only on the lattice grid, the image in the field space in the
present context will just be a set of points in the field space. A continuous
path needs to be drawn in the field space such that the points corresponding
to spins at NN points in physical space are joined by a continuous line. We
assume that the field variation between NN points follows the Geodesic rule
and accordingly draw segments of the trajectory in the field space. We expect
this assumption to be valid for the temperatures we consider as there is
always non-zero correlation between NN points. We mention here that the
Geodesic rule is taken into account while discretising a continuum theory on
the lattice. In numerical simulations of topological defects this rule is used
locate as well as compute the windings of the defects.

In the present case the field space or the OPS is a circle. Note that this
circle can be characterised by an angular variable $\theta$ which can take any
value in $[0-2\pi]$ with the identification of $\theta=0$ and $\theta=2\pi$.
The field configuration along a straight line from ${\rm\bf r}_{i}$ to
${\rm\bf r}_{j}$ is mapped to a trajectory on the circle. The spins
$\vec{s}_{i}$ and $\vec{s}_{j}$ are mapped to the two end points, $\theta_{i}$
and $\theta_{j}$ respectively on the OPS. If the Geodesic rule were valid then
the length, $\Delta\theta_{ij}$, of this trajectory will be minimum of
$|\theta_{i}-\theta_{j}|$ and $2\pi-|\theta_{i}-\theta_{j}|$, i.e
$\Delta\theta_{ij}=\rm{min}[|\theta_{i}-\theta_{j}|,2\pi-|\theta_{i}-\theta_{j}|]$.

The length of the trajectory corresponding to the field configuration from
${\rm\bf r}_{i}$ to ${\rm\bf r}_{j}$ is given by the absolute value of,

$\eta_{ij}(r)=\sum_{k=i}^{j-1}\alpha_{k}\Delta\theta_{k+1,k},~{}\Delta\theta_{k+1,k}=\rm{min}\left[|\theta_{k+1}-\theta_{k}|,2\pi-|\theta_{k+1}-\theta_{k}|\right].$
(5)

$\alpha_{k}=+1(-1)$, if following the Geodesic rule the trajectory traces a
path clockwise(anti-clockwise) in the field space as one goes from the lattice
site ${\rm\bf r}_{k}$ to ${\rm\bf r}_{k+1}$. Note that $\eta_{ij}(r)$
corresponds to net variation of angle $\theta$ of the spins over $r$
separation. In case $\Delta\theta_{ij}\neq\eta_{ij}$, it is counted as the
deviation from the geodesic rule. This computation is repeated for different
values of $r=2,3,.....,\xi_{KM}$ resulting in the function $d(r)$ which gives
the deviation from the Geodesic rule at separation $r$.

In figure 4 we show an example of the field configurations along a linear
stretch of $\xi_{KM}$ for $T=1.43$. $\xi_{KM}=12$ for this $T$. The
corresponding image in the OPS is shown in figure 4. The two end points
correspond to $\theta_{1}\simeq 0.325\pi$ and $\theta_{13}\simeq 0.292\pi$,
with $\Delta\theta_{1,13}=0.033\pi$. $\eta_{1,13}$ in this case turns out to
be $2\pi-\Delta\theta_{1,13}=1.967\pi$; which does not correspond to the
shortest path. This configuration is an example where the Geodesic rule is
violated. In figures 6 another configuration is shown for the same $T$. The
corresponding image map is shown in 6. The endpoints for this configuration
are; $\theta_{1}\simeq 1.76\pi$ and $\theta_{13}\simeq 0.05\pi$. In this case
$\Delta\theta_{1,13}=.29\pi$ but $\eta_{1,13}$ is found to be $2.29\pi$. This
trajectory not only violates the Geodesic rule but also winds the field space,
i.e OPS, once. This implies that for given two end points in the field space
there are many ways in which geodesic rule can be violated but there is only
one way it is followed.

Figure 3: A configuration spins over length $\xi$ in physical space.

Figure 4: The corresponding trajectory $\eta_{1,13}$ defined in Eq.(5) on the
OPS.

Figure 5: A configuration spins over length $\xi$ in physical space.

Figure 6: The corresponding trajectory $\eta_{1,13}$ defined in Eq.(5) on the
OPS.

The Fig.8 and Fig.8 show the probability of deviation from the geodesic rule
vs $r$ at various temperatures in two and three dimensional systems
respectively. For simplicity this calculation was carried out along the $x$
and $y$ directions only. The results were found to be independent of
directions. The results show that the probability of deviation increases with
separation $r$ and temperature. The upper points on the curves for different
temperatures correspond to the deviation at $r=\xi_{KM}$. Note that the
deviation at $r=\xi_{KM}$ increases as we approach $T_{c}$ from above. This is
because $\xi_{KM}$ diverges at $T_{c}$. The larger the $\xi_{KM}$, there are
more possibilities to trace different paths between two end points. It is
observed that paths with multiple windings around the OPS are possible. For
comparison the deviation at the separations $r=\xi$ is also shown in these
figures, indicated by blue dots on the different curves. Though, this
deviation is smaller than that at $r=\xi_{KM}$ but is still significant.

Figure 7: Geodesic rule devation in $2D-$XY model with lattice distance.

Figure 8: Geodesic rule deviation in $3D-$XY model with lattice distance.

###  II.3  Distribution of vortices in two spatial dimensions

In this section the results on the distribution of vortices are presented for
two spatial dimensions. Each of the configuration in the thermal ensemble,
generated using the Monte Carlo simulations, is analysed to compute the number
of defects and net winding number inside a $100\times 100$ sub-lattice, which
is quarter of the full system. Within this sub-lattice each elementary squares
is scanned for vortices. To determine if there is a defect inside a square the
Geodesic rule is applied to specify field variations between NN points. This
uniquely determines the winding number of the field along the perimeter. The
variations of $\theta$ along a closed path on the lattice is a topological
number, i.e $2\pi n$ for integer n. If the $\theta$ of the spins varies by
$(-)2\pi$ as the perimeter is traced once in the clockwise direction, the
winding is considered to be $(-)1$. The net winding on the full lattice is
zero as the the lattice is effectively a torus ($T^{2}$) due to periodic
boundary conditions imposed along the $x,y$ directions. Fig.10 shows, $N_{D}$,
the average number of vortices and anti-vortices in the whole lattice. For
comparison we compute defects inside each elementary square in three
dimensional lattice ignoring the sign of the winding number. The results are
shown in Fig.10. In both cases the number of vortices decreases with
temperature. In Fig.10 and Fig.10, the estimates using the KM are also
included, which underestimate $N_{D}$ by large margins for temperatures
considered here.

Figure 9: $N_{D}$ from Monte Carlo simulations for $2D-$XY model. The lower
curve corresponds to estimate from KM.

Figure 10: $N_{D}$ from Monte Carlo simulations for $3D-$XY model. The lower
curve corresponds to estimate from KM.

The number of vortices for all temperatures $T>T_{c}$, is way above what is
expected from the Kibble Mechanism. It can be easily shown that the
probability of having a vortex or anti-vortex inside an equilateral triangle
of size $\xi_{KM}$ is $1/4$ Bowick:1992rz . If the Geodesic rule is violated
by $50\%$ then this probability that there is a defect inside the triangle
goes up to $5/8$. This suggests that at finite temperature most of the defects
are created via thermal fluctuations. For example, at $T=1.18$, using the
Geodesic rule or it’s violation by $50\%$, the maximum number of defects
expected inside the sub-lattice is about $\sim 27$ when $\xi_{KM}$ is used for
the estimate. The average number of defects in this case is $\sim 700$. This
discrepancy reduces when standard correlation length, $\xi$, and the
probability that the Geodesic rule is violated at this temperature are
considered.

###  II.4 Correlation between vortices and anti-vortices

The Kibble Mechanism predicts strong correlation among the vortices
Digal:1998ak . It is expected that, given a vortex the probability that the
nearest defect is an anti-vortex is higher than a vortex. This can be deduced
from the Kibble Mechanism in the following way Digal:1998ak ;
Vachaspati:1984dz . The number of defects within an area $A$ is $N\sim
A/\xi^{2}$. On the other hand the perimeter of this area is $\sim A^{1/2}$.
This can be traversed in $N^{1/2}\sim A^{1/2}/\xi$ steps of length $\xi$.
Since for each step the average variation of $\theta$ is $\pi/2$, the net
variation of $\theta$ along the perimeter can be computed using method of
random walk of $N^{1/2}$ steps. Hence it is expected that the standard
deviation is $\sim N^{\nu}$ with $\nu=1/4$. For a completely random
distribution of vortices $\nu=1/2$. In Fig.12 we plot the distribution of net
winding number in the sub-lattice for various temperatures.

Figure 11: The distributions of the net winding number in the sub-lattice at
various temperatures.

Figure 12: The exponent $\nu$ corresponding to the standard deviation in the
net winding number distributions.

The distributions shown in Fig.12 can be fitted well with gaussians. The
exponent $\nu$ is given by the standard deviation of the gaussians. The
resulting values of $\nu$ extracted for various temperatures is shown in
Fig.12. The exponent $\nu$ for $T=1.18$ turns out to be $0.18$. For $T=2.5$,
$\nu\simeq 0.21$. Note that a constant factor is to be taken into account
while comparing these results with the the Kibble Mechanism Digal:1998ak . The
resulting $\nu$ approaches the value $1/4$ for $T\to\infty$. This makes sense
as the geodesic rule is imposed at the scale of lattice spacing and the NN
spins are uncorrelated in this limit. In the following we study the Geodesic
rule in $\phi^{4}$ theory. The results are qualitatively similar to the
present case of $O(2)-$spins.

##  III Geodesic rule in $\phi^{4}$ theory

The lattice action for the $\phi^{4}$ theory with $U(1)$ symmetry is taken to
be,

$S_{\phi}=-\kappa\sum_{i,\mu}\left(\phi^{\dagger}_{i}\phi_{i+\hat{\mu}}+h.c\right)+\sum_{i}\left[{1\over
2}\phi^{\dagger}_{i}\phi_{i}+\lambda\left(\phi^{\dagger}_{i}\phi_{i}-1\right)^{2}\right]$
(6)

where $\mu=x,y,z$. $\hat{\mu}$ represent the unit vector in the $\mu-$th
direction. $\lambda$ represents the coupling for the quartic self-interaction
of the $\phi$ field. $\lambda$ is fixed to $1$ for the simulations. The NN
coupling parameter $\kappa$ plays the role of inverse temperature in this
case. The corresponding partition function is given by

${\cal{Z}}=\int\prod_{i}d\phi_{i}^{*}d\phi_{i}~{}{\rm Exp}[-S_{\phi}].$ (7)

The field $\phi$ being complex can be thought of as a two component vector,
whose magnitude can fluctuate unlike in the case of $O(2)-$spins. In the
thermodynamic limit the system undergoes a $2nd$ order phase transition at
$\kappa_{c}$. The paramagnetic phase persists upto $\kappa_{c}$. For
$\kappa>\kappa_{c}$ the system is found to be in the ferromagnetic phase.

The equilibrium configurations in this model are generated by using the pseudo
heat-bath algorithm. In this algorithm to update $\phi_{i}$ an action $S_{i}$
is considered, which is obtained from $S_{\phi}$ by dropping term which does
not depend on $\phi_{i}$. The resulting action $S_{i}$ is rewritten as a sum
of gaussian in $\phi_{i}$ and a quartic term. The gaussian term is used to
generate a new $\phi_{i}$ which is accepted with probability determined by the
change in the quartic term. This process is repeated to update $\phi$ at all
the lattice sites. To reduce the auto correlation over relaxation method is
used. In addition every $5$th configuration is used for analysis. For details
see reference Bunk:1994xs . The simulations are carried out for $60\times
60\times 60$ $3D$ lattice. The Fig.14 and Fig.14 show two dimensional section
of typical field configurations at $\kappa=0.15$ and $\kappa=0.4$
respectively.

Figure 13: A sampled spin configuration
at $T=2T_{c}$.

Figure 14: A sampled spin configuration
at $T=0.76T_{c}$.

At high temperatures, $\kappa<\kappa_{c}$, the volume and thermal average of
$\phi$ is zero, however $\phi_{i}=0$ is suppressed due to vanishing measure.
It is found that the image trajectory, of any field configuration over a
linear stretch $r$, rarely passes through the origin ($\phi=0$) in the field
space. Therefore the phase $\theta$ of the field $\phi$ can be used for the
computation of the probability of deviation from the Geodesic rule following
the procedures adopted in the case of $O(2)-$spins. $\Delta\theta_{ij}$ and
$\eta_{ij}$ are found as function of the separation $r=|{\bf r}_{i}-{\bf
r}_{j}|$. As in the case of $O(2)$ spins whenever
$\Delta\theta_{ij}\neq\eta_{ij}$ the Geodesic rule is considered to be
violated. The Fig.15 shows the probability of deviation from the Geodesic rule
as function of $r$ for different $\kappa$. The deviation increases with $r$
and temperature $1/\kappa$ similar to the case of $O(2)-$spins. The red
squares on the curves correspond to the deviation at $r=\xi_{KM}$.

Figure 15: Geodesic rule deviation for $\phi^{4}$ theory.

##  IV Conclusions

We have studied the variations of field in the $XY-$model in two and three
dimensions and $\phi^{4}$ theory in three dimensions using Monte Carlo lattice
simulations. For two dimensions the lattice size was taken to be $200\times
200$. For three dimensions we consider $60^{3}$ lattice. The computations are
carried out above the critical temperature to check the validity of the
Geodesic rule used in the Kibble-Zurek mechanism of defect formation in 2nd
order phase transitions. For this purpose a correlation length $\xi_{KM}$ is
defined using the criteria that field correlation vanish approximately. This
length scale is different from the conventional correlation length $\xi$. At
$\xi$ separations the field can not be considered uncorrelated which is
crucial for the Geodesic rule. It is found that even at $\xi$ separations the
field variation deviates significantly from what is expected from the Geodesic
rule. At $\xi_{KM}$ separations non geodesic paths are found to be equally or
more probable than the geodesic paths.

We also studied the number of vortices and the net winding number in the case
of $O(2)-$spins in two dimensions. The geodesic rule is assumed for field
variations between nearest neighbour spins. This is necessary to probe the
location and winding charge of the vortices. The analysis is carried out a
sub-lattice of size $100\times 100$ which is a quarter of the whole lattice.
This is necessary as in this case the net winding number is not constrained to
be zero. The number of vortices is found to be significantly higher than what
is expected from the Kibble-Zurek mechanism. The discrepancy does not improve
much even after the deviation from the geodesic rule is incorporated in the
Kibble-Zurek mechanism. The simulation results show that the large number of
vortices are a result of the thermal production. This is supported by results
of the distribution of the net winding number. The results show that the
corresponding standard deviation is much smaller than what is expected from
the consideration of the Geodesic rule. This result suggests that there is
significant pairing between the vortices and anti-vortices. It is expected
that from topological considerations thermal fluctuations always form pair of
vortex-anti-vortex and at smaller separations. We mention here that in
theories with gauge symmetries there is no satisfactory argument for the
Geodesic rule Rudaz:1992wy . The trajectory of the field between two domains
can be deformed by gauge transformations. The Monte Carlo simulation method
presented here can be extended to settle the Geodesic rule in these theories,
which we plan to do next.

###### Acknowledgements.

We thank A. P. Balachandran, Ajit. M. Srivastava and Rajarshi Ray for valuable
comments.

REFERENCES

## References

  * (1) I. Chuang, B. Yurke, R. Durrer and N. Turok, Science 251, 1336-1342 (1991) doi:10.1126/science.251.4999.1336
  * (2) P.C. Hendry et al, Nature (London) 368, 315 (1994) doi:10.1038/368315a0
  * (3) V.M.H. Ruutu et al, Nature(London) 382, 334 (1996) doi:10.1038/382334a0
  * (4) C. Bäuerle el al., Nature 382 (1996) 332 doi:10.1038/382332a0
  * (5) M. E. Dodd, P. C. Hendry, N. S. Lawson, P. V. E. McClintock and C. D. H. Williams, Phys. Rev. Lett. 81, no.17, 3703-3706 (1998) doi:10.1103/PhysRevLett.81.3703 [arXiv:cond-mat/9808117 [cond-mat.soft]].
  * (6) S. Digal, R. Ray and A. M. Srivastava, Phys. Rev. Lett. 83, 5030 (1999) doi:10.1103/PhysRevLett.83.5030 [arXiv:hep-ph/9805502 [hep-ph]].
  * (7) R. Carmi, E. Polturak and G. Koren, Phys. Rev. Lett. 84, 4966-4969 (2000) doi:10.1103/PhysRevLett.84.4966
  * (8) M. B. Hindmarsh and T. W. B. Kibble, Rept. Prog. Phys. 58, 477-562 (1995) doi:10.1088/0034-4885/58/5/001 [arXiv:hep-ph/9411342 [hep-ph]].
  * (9) J. Magueijo and R. H. Brandenberger, [arXiv:astro-ph/0002030 [astro-ph]].
  * (10) A. M. Srivastava, P. Bagchi, A. Das and B. Layek, Pramana 89, no.4, 68 (2017) doi:10.1007/s12043-017-1465-1
  * (11) S. Digal, R. Ray, S. Sengupta and A. M. Srivastava, Phys. Rev. Lett. 84, 826-829 (2000) doi:10.1103/PhysRevLett.84.826 [arXiv:hep-ph/9911446 [hep-ph]].
  * (12) W. H. Zurek, W. H. Zurek, Nature (London) 317, 505 (1985) doi:10.1038/317505a0; Phys. Rept. 276, 177-221 (1996) doi:10.1016/S0370-1573(96)00009-9 [arXiv:cond-mat/9607135 [cond-mat]].
  * (13) G. Toulouse, M. Kléman. Journal de Physique Lettres, Edp sciences, 1976, 37 (6), pp.149-151. 10.1051/jphyslet:01976003706014900. jpa-00231261
  * (14) N. D. Mermin, Rev. Mod. Phys. 51, 591-648 (1979) doi:10.1103/RevModPhys.51.591
  * (15) Srivastava A.M. (2001) Topological Defects in Condensed Matter Physics. In: Field Theories in Condensed Matter Physics. Texts and Readings in Physical Sciences doi:10.1007/978-93-86279-07-1 5
  * (16) Paul M Chaikin, T C Lubensky; Principles of condensed matter physics, Cambridge,Cambridge University Press,1995, ISBN:0 521 43224 3
  * (17) T. W. B. Kibble, J. Phys. A 9, 1387-1398 (1976) doi:10.1088/0305-4470/9/8/029
  * (18) T. W. B. Kibble and A. Vilenkin, Phys. Rev. D 52, 679-688 (1995) doi:10.1103/PhysRevD.52.679 [arXiv:hep-ph/9501266 [hep-ph]].
  * (19) M. J. Bowick, L. Chandar, E. A. Schiff and A. M. Srivastava, Science 263, 943-945 (1994) doi:10.1126/science.263.5149.943 [arXiv:hep-ph/9208233 [hep-ph]].
  * (20) R. Ray and A. M. Srivastava, Phys. Rev. D 69, 103525 (2004) doi:10.1103/PhysRevD.69.103525 [arXiv:hep-ph/0110165 [hep-ph]].
  * (21) A. M. Srivastava, Phys. Rev. D 46, 1353-1367 (1992) doi:10.1103/PhysRevD.46.1353
  * (22) S. Chakravarty and A. M. Srivastava, Nucl. Phys. B 406, 795-807 (1993) doi:10.1016/0550-3213(93)90010-M [arXiv:hep-ph/9209246 [hep-ph]].
  * (23) E. J. Copeland and P. M. Saffin, Phys. Rev. D 54, 6088-6094 (1996) doi:10.1103/PhysRevD.54.6088 [arXiv:hep-ph/9604231 [hep-ph]].
  * (24) S. Digal, S. Sengupta and A. M. Srivastava, Phys. Rev. D 56, 2035-2043 (1997) doi:10.1103/PhysRevD.56.2035 [arXiv:hep-ph/9705246 [hep-ph]].
  * (25) S. Digal and A. M. Srivastava, Phys. Rev. Lett. 76, 583-586 (1996) doi:10.1103/PhysRevLett.76.583 [arXiv:hep-ph/9509263 [hep-ph]].
  * (26) J. M. Kosterlitz and D. J. Thouless, J. Phys. C 6, 1181-1203 (1973) doi:10.1088/0022-3719/6/7/010
  * (27) Hasting W.K., “Monte Carlo sampling methods using Markov chains and their applications”; Biometrika 57 (1970) 97–109, doi:10.1093/biomet/57.1.97
  * (28) T. Vachaspati and A. Vilenkin, Phys. Rev. D 30, 2036 (1984) doi:10.1103/PhysRevD.30.2036
  * (29) B. Bunk, Nucl. Phys. B Proc. Suppl. 42, 566-568 (1995) doi:10.1016/0920-5632(95)00313-X
  * (30) S. Rudaz and A. Mohan Srivastava, Mod. Phys. Lett. A 8, 1443-1450 (1993) doi:10.1142/S0217732393001161 [arXiv:hep-ph/9212279 [hep-ph]].


# Electrochemical oxidation of CO on Cu single crystals under alkaline
conditions

Aarti Tiwari  Department of Physics, Technical University of Denmark (DTU),
Fysikvej 311, 2800 Kgs. Lyngby, Denmark  Hendrik H. Heenen  Department of
Physics, Technical University of Denmark (DTU), Fysikvej 311, 2800 Kgs.
Lyngby, Denmark  Anton Simon Bjørnlund  Department of Physics, Technical
University of Denmark (DTU), Fysikvej 311, 2800 Kgs. Lyngby, Denmark
Degenhart Hochfilzer  Department of Physics, Technical University of Denmark
(DTU), Fysikvej 311, 2800 Kgs. Lyngby, Denmark  Karen Chan  Department of
Physics, Technical University of Denmark (DTU), Fysikvej 311, 2800 Kgs.
Lyngby, Denmark<EMAIL_ADDRESS>Sebastian Horch  Department of Physics,
Technical University of Denmark (DTU), Fysikvej 311, 2800 Kgs. Lyngby, Denmark
<EMAIL_ADDRESS>
###### Abstract

We perform a joint experimental-theoretical study of the electrochemical
oxidation of CO on copper (Cu) under alkaline conditions. Using cyclic
voltammetry on Cu single crystal surfaces, we demonstrate that both Cu
terraces and steps show CO oxidation activity at potentials just slightly
positive (0.03–0.14 V) of the thermodynamic equilibrium potential. The
overpotentials are 0.23–0.12 V lower than that of gold ($\approx 0.26$ V),
which up until now has been considered to be the most active catalyst for this
process. Our theoretical calculations suggest that Cu’s activity arises from
the advantageous combination of simultaneous *OH adsorption under CO oxidation
potentials and surmountable *CO-*OH coupling barriers. Experimentally observed
onset potentials are in agreement with the computed onsets of *OH adsorption.
We furthermore show that the onsets of *OH adsorption on steps are more
affected by *CO-*OH interactions than on terraces due to a stronger
competitive adsorption. Overall, Cu(100) shows the lowest overpotential (0.03
V) of the facets considered.

###### keywords:

American Chemical Society, LaTeX

These authors contributed equally to this work These authors contributed
equally to this work These authors contributed equally to this work These
authors contributed equally to this work These authors contributed equally to
this work Center for Bachelor of Engineering Studies, Technical University of
Denmark (DTU), Lautrupvang 15, 2750 Ballerup, Denmark Cu, CV, SC

In this paper, we revisit the electrochemical (EC) interaction of copper (Cu)
with CO at oxidative potentials under alkaline conditions. The expected
reaction, the EC oxidation of CO on Cu seems to have been overlooked since the
early works of Kita _et al._ in the 1980s, probably because they concluded
that Cu was inactive for CO oxidation 1. This conclusion resulted in most
subsequent research being focused on Au, leading to the understanding that it
is the most active catalyst for EC CO oxidation 1, 2, 3, 4, 5. When we tried
to reproduce Kita’s measurements, we did, however, find oxidative processes on
Cu induced by the presence of CO already at potentials close to the
thermodynamic potential for CO oxidation. This result suggests Cu as a very
good catalyst for the electrochemical oxidation of CO. Cu could thus be an
interesting, cheap alternative to the known expensive catalysts for EC CO
oxidation e.g. in fuel cells, where CO has to be eliminated from the feed to
avoid poisoning of especially Platinum catalysts.6

In a joint experimental-theoretical effort, we explore the interaction of Cu
with *CO and *OH for alkaline CO oxidation. To this end, we have measured
cyclic voltammograms (CVs) with and without CO on well-defined Cu single
crystals (SCs). These CVs reveal oxidative processes at the Cu/electrolyte
interface in the presence of CO which we interpret as the EC oxidation of CO.
We corroborate this hypothesis via _ab initio_ simulations demonstrating a
feasible reaction mechanism and interface composition at these potentials. We
studied Cu SCs with four different surface structures, namely the low-index
Cu(100), (110), (111) and (211) orientations. The EC measurements were
performed in the respective Cu SC’s fingerprint potential range; this
potential range is positive of the onset of hydrogen evolution and negative of
the surface oxidation region. The fingerprint region therefore corresponds to
a potential range where the Cu SC is not an oxide or (potentially) a hydride7.
The characteristic fingerprint EC response of the first three orientations,
which expose large terraces and relatively few steps in absence of any
reactive species are known from a previous study (details in Methods section)
8. Here, we for the first time also present the fingerprint response of (211)
(see Figure S1, SI), a prototype of a stepped surface (as further detailed in
the SI Sec. II).

Figure 1 shows CVs of the four Cu orientations in both blank (Ar-) and CO-
saturated electrolyte. These measurements show that (a) all four facets have a
distinct fingerprint response in blank electrolyte (grey traces), and (b) the
fingerprint CVs show clear oxidative features due the presence of CO in the
electrolyte (colored traces). In the presence of CO, all facets show increased
anodic currents at potentials positive of $-0.1$ V (vs. RHE as described
further down in Methods section) which indicates an oxidation process (_onset_
indicated by cyan dotted lines). These currents appear slightly positive of
both OH-adsorption (known from the fingerprint CVs,8 _reversible potential_
indicated by vertical black dotted lines for better visibility) and the
standard equilibrium potential for the CO oxidation reaction
$(\ch{CO}+\ch{H2O}\rightarrow\ch{CO2}+\ch{2H^{+}}+\ch{2e^{-}})$ at $-0.10$ V
1. In Table 1, we list the experimental onset potentials for the reaction,
which we define as the potential corresponding to a current density of 25
$\mu$A cm-2. The oxidative features vanish as soon as CO is purged out of the
electrolyte and after which the fingerprint CVs corresponding to an Ar-purged
environment are recovered. This reversibility indicates that the surfaces have
not undergone any irreversible change in the presence of CO.

Figure 1: Cyclic voltammograms for Cu(211), (110), (100) and (111) in the
respective fingerprint regions measured in blank Ar-saturated (grey trace) and
CO-saturated (colored trace) 0.1 M KOH at a scan rate of 50 mV/s. Vertical
black dotted lines represent the reversible OH adsorption potential in absence
of CO whereas cyan dotted lines mark the onset of CO oxidation.

To investigate the reaction mechanism and the interplay of *CO and *OH
adsorption, we performed _ab initio_ simulations based on density functional
theory (DFT) of the EC response of the four Cu facets to these adsorbates at
anodic potentials. We considered a Langmuir-Hinshelwood mechanism as follows:

$\displaystyle\mathrm{CO(g)+\,^{*}}$
$\displaystyle\rightleftharpoons\mathrm{{}^{*}CO}$ (1)
$\displaystyle\mathrm{H_{2}O(l)+\,^{*}}$
$\displaystyle\rightleftharpoons\mathrm{{}^{*}OH+H^{+}+e^{-}}$ (2)
$\displaystyle\mathrm{{}^{*}CO+\,^{*}OH}$
$\displaystyle\rightleftharpoons\mathrm{{}^{*}COOH}$ (3)
$\displaystyle\mathrm{{}^{*}COOH+\,^{*}OH}$
$\displaystyle\rightleftharpoons\mathrm{H_{2}O(l)+CO_{2}(g)+2^{*}}$ (4)

where * denotes an empty site on the Cu surface. The corresponding free energy
diagrams are shown in Figure 2, and we have explicitly considered the
association barriers for steps 3 and 4. In this mechanism, *OH adsorption (2)
is the sole electrochemical step and the chemical association of *CO and *OH
(3) is the rate-determining step on all facets considered. We show in Figure 2
the energetics both at the theoretically calculated equilibrium potential
(-0.22V${}_{\text{RHE}}$, dashed lines), as well as at a theoretical onset
potential ($E\mathrm{{}^{Th}_{onset}}$ solid lines), which we define to
correspond to the potential where *OH adsorption becomes exergonic. At
$E\mathrm{{}^{Th}_{onset}}$, almost all potential-independent reaction
thermodynamics and barriers are $<0.75$ eV and therefore easily surmountable
at room temperature 9 (c.f. Figure 2). Only on the step-edge of (211), we find
the *CO-*OH coupling barrier to be 0.75 eV and therefore a room-temperature
activity at the step-edge is equivocal. On the corresponding single crystal,
CO oxidation would nevertheless still be feasible on its (111) terrace which
could also give rise to the observed current (c.f. SI Sec. IV). In the
proposed reaction mechanism, the Faradaic current arises through continuous
replenishment of *OH on all four Cu surfaces which is consumed in the CO
oxidation reaction.

Figure 2: Free energy diagram for the suggested CO oxidation reaction
mechanism on Cu(111), Cu(100), Cu(110) and Cu(211). The free energy profile is
depicted for the onset potential (solid lines) and for the thermodynamic
equilibrium potential (dashed lines) as predicted by DFT. The free energies of
the chemical (i.e. potential-independent) barriers are annotated in eV.

In our suggested mechanism, the electrochemical response only depends on *OH
adsorption, since this is the only involved electrochemical step. Thus, the
facet-dependent adsorption of *OH determines the onset of the EC CO oxidation.
Indeed, as shown in Table 1, the theoretically calculated onset potentials –
solely depending on the interaction-free *OH adsorption energies – agree
within DFT-accuracy with the experimental values. The larger deviation between
the theoretical and experimental onset potentials for the large-terrace facets
may originate in an underestimated solvation stabilization of *OH through the
implicit solvent, as recently demonstrated in Ref. 10 . As shown in Table 1,
applying an improved solvation correction (c.f. SI Sec. III) leads to a
considerably better agreement.

The experimental onset potentials suggests (100) to be the most active facet
towards CO oxidation with an overpotential of 0.03 V which is ca. 0.23 V less
than Au, previously thought to be the most electroactive element for CO
oxidation. 1 Our proposed mechanism rationalizes this exceptionally low
overpotential of CO oxidation on Cu. The potentials at which *OH is predicted
to adsorb on the Cu facets and thus CO oxidation sets on (c.f. Figure 2) are
very close to the thermodynamic equilibrium potential of CO oxidation
(calculated to be -0.22 V which matches the experimentally determined value of
-0.1 V within DFT accuracy). This oxophilicity of Cu is in contrast to Au or
Pt, for example, where *OH adsorption occurs significantly more positive 10 of
the thermodynamic equilibrium potential of CO oxidation, thus leading to
higher overpotentials. Due to very high adsorption potentials of *OH in the
case of Au, an alternative Eley-Rideal mechanism involving the coupling of *CO
with OH- in solution has been suggested 3, 4, 11.

Table 1:  Experimental onset potentials ($E\mathrm{{}_{onset}}$), defined to be at 25 $\mu$A cm-2, as well as the theoretical onset potentials, defined as the onset of *OH adsorption, without ($E\mathrm{{}_{onset}^{Th}}$) and with ($\tilde{E}\mathrm{{}_{onset}^{Th}}$) solvation corrections (c.f. Fig. 2). |  $E\mathrm{{}_{onset}}$ (V) |  $E\mathrm{{}_{onset}^{Th}}$ (V) |  $\tilde{E}\mathrm{{}_{onset}^{Th}}$ (V)
---|---|---|---
Cu(111) | 0.04 | 0.36 | 0.18
Cu(100) | -0.07 | 0.16 | 0.14
Cu(110) | 0.02 | -0.04 | 0.02
Cu(211) | -0.05 | -0.11 | -0.03

In general, Figure 1 shows that, on the facets exposing large terraces (100
and 111) the onsets of CO oxidation (cyan dotted lines) and reversible *OH
adsorption (black dotted lines) in absence of CO essentially coincide to
within 0.1 V. In contrast, on the step rich facets (211 and 110) CO oxidation
occurs positive of the onset of *OH adsorption on the step-edges in an Ar
atmosphere which we find at -0.35 V for (110) (c.f. Figure 1) and
theoretically estimate to be $\leq$ -0.4 V for (211) (see SI Sec. II).

We attribute the positive potential shift between CO oxidation and *OH
adsorption in the absence of CO to arise from the stronger competitive
adsorption of CO vs. OH on the steps. This effect is apparent from the
theoretical surface Pourbaix diagrams, shown in Figure 3. From the
configurations sampled, the onset of *OH adsorption in the presence of *CO is
shifted +0.14 V on the terraces, and +0.21/+0.42 V on the (110)/(211) steps.
This competitive effect can also be seen through the additional/shifted peak
at -0.15 V in a CO atmosphere on 110. Following our DFT simulations we,
however, do not assign this feature to the *OH adsorption on the 110 step-edge
but to a different site, e.g. on the 1$\times$2 reconstructed (110) facet 8.
Figure 1 does not depict an *OH adsorption feature for the step-edge on (211)
under Ar atmosphere. Our simulations suggest *OH adsorption on the (211) step
edge to occur $\leq$ -0.4 V vs. RHE, and therefore the adsorption peak would
be convoluted with HER currents (c.f. SI Sec 2). Nevertheless, without the
presence of CO, we still observe a broad 0.12 V peak which corresponds to *OH
adsorption on the corresponding (111) terrace of the (211) single crystal (see
Fig. S1, SI) 8. In presence of CO, no *OH adsorption features can be
identified on (211). We believe these features, which are subjected to
adsorption competition in presence of CO and would thus be found at more
positive potentials, to be hidden within the CO oxidation onset.

We note two observations on the potential dependence of the Faradaic currents
(c.f. Figure 1):

  1. 1.

On (111), (100) and (110), the Faradaic current is almost constant with
respect to potential. This can be rationalized by the presence of a potential
independent rate determining step and relatively constant coverages of the
adsorbates involved. Alternatively, a CO diffusion limitation can also lead to
a constant current profile.

  2. 2.

On (211) such a potential independent current is only seen from -0.01 to 0.1 V
and a second current with an apparent potential dependence sets in at $\geq
0.1$ V that also levels off as potential increases. We attribute the first
Faradaic process to oxidation of *CO adsorbates on the (111) terrace portion
on (211), where the chemical *CO-*OH coupling is more feasible. The onset of
this current is more positive than for the (111) single crystal, which
suggests that the *OH is initially bound to the stronger binding (211) step
edge (see SI Sec. IV for a rationalized pathway). We speculate that the second
Faradaic process may also arise from the oxidation of CO on (211), but at a
different active site (i.e. the step edge where the direct *OH/*CO coupling is
not very facile). Due to its apparent potential dependence, this alternative
process may arise from an electrochemical rate determining step (i.e. *CO +
OH- $\rightarrow$ *COOH + $e^{-}$).

Figure 3: Pourbaix analysis of the adsorbates H* (green), *OH (red), *CO
(blue), O* (purple), and co-adsorbed *CO/*OH (orange) on Cu(111), Cu(100),
Cu(110), and Cu(211). The black solid line represents the clean Cu surface.
The dashed lines represent the most stable configuration of each adsorbate
(over varying coverage) and the labeled solid lines the overall most stable
configuration over the potential range E. The energy $\Delta
G_{ads}^{N_{\theta}}$ is the adsorption energy normalized to the coverage
(i.e. per site). The shaded area corresponds to the potential range in which
the coverage of *CO and *OH are $\geq$ 5% (dark grey) and $\geq$ 1% (light
grey) as estimated from a Boltzmann distribution (see text). The dotted lines
indicate the theoretical interaction-free shift in *OH adsorption potential
through the presence of *CO.

The proposed mechanism necessitates a finite coverage of both *CO and *OH on
the Cu surfaces for the chemical reaction steps to take place. This condition
must therefore be fulfilled at all potentials where CO oxidation is active
(i.e. $\approx 0.0\ldots 0.4$ V). We investigate the composition of surface
coverage by a Pourbaix analysis 12, 13. For this analysis, we sample the
adsorbates O*, H*, *OH, and *CO on the individual Cu facets. To introduce a
coverage dependence, we include each adsorbate in a series of different
supercells and further conduct a configurational sampling for *OH and *CO and
mixed *OH/*CO (which are dominating at the potentials of interest) in confined
supercells enabling (mixed) coverages from 25–100 % (see details SI Sec. III).
We approximate the potential dependent adsorbate coverages via the sampled
configurations for each surface using a Boltzmann distribution. From the most
stable free energies $G_{i}(E)$ of each surface configuration $i$ with an
adsorbate coverage $\theta_{i}^{ads}$, we compute the total coverage
$\theta_{tot}^{ads}$ at the potential $E$ for each adsorbate $ads$:

$\displaystyle\theta_{tot}^{ads}(E)=\sum_{i}^{N}\theta_{i}^{ads}\exp\left(\frac{G_{i}(E)}{\mathrm{k_{B}}T}\right)/Z(E)$
$\displaystyle\mathrm{with}\;Z(E)=\sum_{i}^{N}\exp\left(\frac{G_{i}(E)}{\mathrm{k_{B}}T}\right)$

where $N$ is the number of sampled compositions and $Z(E)$ the potential
dependent partition function. While providing indicative potential dependent
surface coverages, we note that the computed partition function is not
exhaustive since our sampling is conducted in a limited configurational space.
Additionally, the specific interaction with water may affect the considered
adsorbate-adsorbate interactions 14 and could lead to slight deviation from an
idealized behaviour assumed in our analysis. The shaded areas in Figure 3
highlight the potentials at which the coverages of *CO
($\theta_{i}^{*\mathrm{CO}}$) and *OH ($\theta_{i}^{*\mathrm{OH}}$) are $\geq$
5% or $\geq$ 1%, respectively, according to our Pourbaix analysis. It can be
clearly seen, that finite coverages of both *CO and *OH are present at the
potentials where CO oxidation occurs. We note that, while *OH and *CO compete
for adsorption sites – as evident by the shifted *OH adsorption peaks in the
CV spectra (c.f. Figure 1) – we nevertheless find them to moderately co-
stabilize in mixed coverage regimes.

In regard of the suggested Langmuir-Hinshelwood mechanism, the following
observations can finally be summarized in its support:

  1. 1.

coexistence of *OH and *CO at the relevant potential ranges

  2. 2.

low *OH-*CO coupling barriers

  3. 3.

the CO oxidation onset potential coincides with the *OH-adsorption (as
theoretically predicted and experimentally confirmed)

Further, a Langmuir-Hinshelwood mechanism and a chemical reaction as a rate
determining step is also in line with previous reports of CO oxidation on
Pt(111) 15. We, nevertheless, note that an alternative Eley-Rideal mechanism
involving a direct oxidation of adsorbed *CO through H2O cannot be definitely
excluded.

In summary, we experimentally and theoretically show that Cu is active for the
EC oxidation of CO. We find that the corresponding overpotentials are
surprisingly low compared to Au, which has been considered as the most
electroactive catalyst so far under alkaline conditions. Furthermore, our
theoretical calculations suggest that Cu’s superior overpotential can be
attributed to the facile adsorption of *OH on the catalyst surface as well as
the surmountable *CO-*OH coupling barriers. We suggest that Cu is more active
than Au or Pt because the latter adsorb *OH at much more positive potentials.
Additionally, a structural dependence towards CO oxidation was observed
between the four measured orientations, which arises from the facet
specificity of *OH and *CO co-adsorption. Our findings highlight the
electrocatalytic capability of Cu for CO oxidation due to its fine balance of
adsorption energies for the involved adsorbates. Further, this work suggests
Cu as a promising low-cost candidate for EC purification of CO from gas feed
in fuel cells.

Methods

Measurements were performed in a conventional EC setup (described elsewhere
16, 8) in either Ar- or CO-saturated 0.1 M KOH electrolyte and all potentials
are referred against reversible hydrogen electrode (RHE). The SCs were cleaned
by electropolishing 8 in a 66% H3PO4 (85% EMSURE, Merck) at 2.0 V for 60 s
followed by rinsing with Millipore water (18.2 M$\Omega$cm) prior to recording
the CVs. The cleanliness of the SCs was confirmed on a regular basis with _ex
situ_ XPS after the measurements and no traces of contamination especially
metals were found. The potential ranges for the CV measurements for each
fingerprint region in absence of CO in the electrolyte were purposely chosen
in order to avoid the extremes of hydrogen evolution and surface oxidation to
allow focus on the structure and potential dependence of the characteristic OH
adsorption/desorption features.8 All the measurements have been performed at
least thrice and on two independent sets of Cu SCs to ensure reproducibility.

Density functional theory (DFT) based on the BEEF-vdW 17 exchange-correlation
functional as implemented in Quantum-Espresso 18 was used. The electrochemical
environment was mimicked through an implicit solvent 19, 20 and
electrochemical reactions were referenced via the computational hydrogen-
electrode 21. DFT energies were corrected for CO gas-phase reaction
thermodynamics 22. Atomic structures, geometry optimizations and nudged-
elastic-band (NEB) calculations 23, 24 were handled via ASE25, 26, 27, 28, 29
(see further details in the SI Sec. III).

Note added in proof

We were made aware of a parallel publication on CO oxidation over Cu(111) in
alkaline electrolyte 30 that appeared during the review of this paper. Our
experimental results are consistent with their work and their electrochemical
infrared spectroscopy data supports our mechanistic description. However, we
have different theoretical interpretations of the active site and Cu’s
inherent activity.

Separate representation of the Ar and CO saturated CV spectra; detailed
discussion of *OH adsorption peaks on Cu(211); details of DFT calculations;
detailed depiction of reaction barriers for CO oxidation reactions.

This work was supported by research grant 9455 from VILLUM FONDEN and the
European Union’s Horizon 2020 research and innovation programme under the
Marie Sklodowska-Curie grant agreement no. 713683.

## References

  * Kita et al. 1985 Kita, H.; Nakajima, H.; Hayashi, K. Electrochemical oxidation of CO on Au in alkaline solution. _Journal of Electroanalytical Chemistry and Interfacial Electrochemistry_ 1985, _190_ , 141–156
  * Chang et al. 1991 Chang, S. C.; Hamelin, A.; Weaver, M. J. Dependence of the electrooxidation rates of carbon monoxide at gold on the surface crystallographic orientation: a combined kinetic-surface infrared spectroscopic study. _The Journal of Physical Chemistry_ 1991, _95_ , 5560–5567
  * Edens et al. 1996 Edens, G. J.; Hamelin, A.; Weaver, M. J. Mechanism of Carbon Monoxide Electrooxidation on Monocrystalline Gold Surfaces: Identification of a Hydroxycarbonyl Intermediate. _The Journal of Physical Chemistry_ 1996, _100_ , 2322–2329
  * Blizanac et al. 2004 Blizanac, B. B.; Lucas, C. A.; Gallagher, M. E.; Arenz, M.; Ross, P. N.; Marković, N. M. Anion Adsorption, CO Oxidation, and Oxygen Reduction Reaction on a Au(100) Surface: The pH Effect. _Journal of Physical Chemistry B_ 2004, _108_ , 625–634
  * Rodriguez et al. 2010 Rodriguez, P.; Garcia-Araez, N.; Koverga, A.; Frank, S.; Koper, M. T. CO electroxidation on gold in alkaline media: A combined electrochemical, spectroscopic, and DFT study. _Langmuir_ 2010, _26_ , 12425–12432
  * Valdés-López et al. 2020 Valdés-López, V. F.; Mason, T.; Shearing, P. R.; Brett, D. J. Carbon monoxide poisoning and mitigation strategies for polymer electrolyte membrane fuel cells – A review. _Progress in Energy and Combustion Science_ 2020, _79_ , 100842
  * Huang et al. 2018 Huang, J.; Hörmann, N.; Oveisi, E.; Loiudice, A.; De Gregorio, G. L.; Andreussi, O.; Marzari, N.; Buonsanti, R. Potential-induced nanoclustering of metallic catalysts during electrochemical CO 2 reduction. _Nat. Commun._ 2018, _9_ , 1–9
  * Tiwari et al. 2020 Tiwari, A.; Heenen, H. H.; Bjørnlund, A. S.; Maagaard, T.; Cho, E.; Chorkendorff, I.; Kristoffersen, H. H.; Chan, K.; Horch, S. Fingerprint Voltammograms of Copper Single Crystals under Alkaline Conditions: A Fundamental Mechanistic Analysis. _J. Phys. Chem. Lett._ 2020, _11_ , 1450–1455
  * Nørskov et al. 2014 Nørskov, J. K.; Studt, F.; Abild-Pedersen, F.; Bligaard, T. _Fundamental Concepts in Heterogeneous Catalysis_ ; John Wiley & Sons, Inc., 2014
  * Heenen et al. 2020 Heenen, H. H.; Gauthier, J. A.; Kristoffersen, H. H.; Ludwig, T.; Chan, K. Solvation at metal/water interfaces: An ab initio molecular dynamics benchmark of common computational approaches. _J. Chem. Phys._ 2020, _152_ , 144703
  * Duan and Henkelman 2018 Duan, Z.; Henkelman, G. Calculations of the pH-Dependent Onset Potential for CO Electrooxidation on Au (111). _Langmuir_ 2018, _34_ , 15268–15275
  * Hansen et al. 2008 Hansen, H. A.; Rossmeisl, J.; Nørskov, J. K. Surface Pourbaix diagrams and oxygen reduction activity of Pt, Ag and Ni (111) surfaces studied by DFT. _Phys. Chem. Chem. Phys._ 2008, _10_ , 3722–3730
  * Bukas et al. 2018 Bukas, V. J.; Kim, H. W.; Sengpiel, R.; Knudsen, K.; Voss, J.; McCloskey, B. D.; Luntz, A. C. Combining experiment and theory to unravel the mechanism of two-electron oxygen reduction at a selective and active co-catalyst. _ACS Catal._ 2018, _8_ , 11940–11951
  * Kristoffersen et al. 2018 Kristoffersen, H. H.; Vegge, T.; Hansen, H. A. OH formation and H2 adsorption at the liquid water–Pt(111) interface. _Chem. Sci._ 2018, _9_ , 6912–6921
  * Lebedeva et al. 2002 Lebedeva, N.; Koper, M.; Feliu, J.; Van Santen, R. Mechanism and kinetics of the electrochemical CO adlayer oxidation on Pt (111). _J. Electroanal. Chem._ 2002, _524_ , 242–251
  * Tiwari et al. 2019 Tiwari, A.; Maagaard, T.; Chorkendorff, I.; Horch, S. Effect of Dissolved Glassware on the Structure-Sensitive Part of the Cu(111) Voltammogram in KOH. _ACS Energy Lett._ 2019, _4_ , 1645–1649
  * Wellendorff et al. 2012 Wellendorff, J.; Lundgaard, K. T.; Møgelhøj, A.; Petzold, V.; Landis, D. D.; Nørskov, J. K.; Bligaard, T.; Jacobsen, K. W. Density functionals for surface science: Exchange-correlation model development with Bayesian error estimation. _Phys. Rev. B_ 2012, _85_ , 235149
  * Giannozzi et al. 2009 Giannozzi, P. et al. QUANTUM ESPRESSO: a modular and open-source software project for quantum simulations of materials. _J. Phys. Condens. Matter_ 2009, _21_ , 395502
  * Andreussi et al. 2012 Andreussi, O.; Dabo, I.; Marzari, N. Revised self-consistent continuum solvation in electronic-structure calculations. _J. Chem. Phys._ 2012, _136_
  * 20 Andreussi, O. See http://www.quantum-environment.org for ENVIRON package. http://www.quantum-environment.org
  * Nørskov et al. 2004 Nørskov, J.; Rossmeisl, J.; Logadottir, A.; Lindqvist, L.; Kitchin, J.; Bligaard, T.; Jónsson, H. Origin of the Overpotential for Oxygen Reduction at a Fuel-Cell Cathode. _J. Phys. Chem. B_ 2004, _108_ , 17886–17892
  * Christensen et al. 2015 Christensen, R.; Hansen, H. A.; Vegge, T. Identifying systematic DFT errors in catalytic reactions. _Catal. Sci. Technol._ 2015, _5_ , 4946–4949
  * Henkelman et al. 2000 Henkelman, G.; Uberuaga, B. P.; Jónsson, H. A climbing image nudged elastic band method for finding saddle points and minimum energy paths. _J. Chem. Phys._ 2000, _113_ , 9901–9904
  * Henkelman and Jónsson 2000 Henkelman, G.; Jónsson, H. Improved tangent estimate in the nudged elastic band method for finding minimum energy paths and saddle points. _J. Chem. Phys._ 2000, _113_ , 9978–9985
  * Larsen et al. 2017 Larsen, A. H. et al. The atomic simulation environment—a Python library for working with atoms. _J. Phys. Condens. Matter_ 2017, _29_ , 273002
  * Torres et al. 2019 Torres, J. A. G.; Jennings, P. C.; Hansen, M. H.; Boes, J. R.; Bligaard, T. Low-scaling algorithm for nudged elastic band calculations using a surrogate machine learning model. _Phys. Rev. Lett._ 2019, _122_ , 156001
  * Koistinen et al. 2017 Koistinen, O.-P.; Dagbjartsdóttir, F. B.; Ásgeirsson, V.; Vehtari, A.; Jónsson, H. Nudged elastic band calculations accelerated with Gaussian process regression. _J. Chem. Phys._ 2017, _147_ , 152720
  * del Río et al. 2019 del Río, E. G.; Mortensen, J. J.; Jacobsen, K. W. Local Bayesian optimizer for atomic structures. _Phys. Rev. B_ 2019, _100_ , 104103
  * Boes et al. 2019 Boes, J. R.; Mamun, O.; Winther, K.; Bligaard, T. Graph Theory Approach to High-Throughput Surface Adsorption Structure Generation. _J. Phys. Chem. A_ 2019, _123_ , 2281–2285
  * Auer et al. 2020 Auer, A.; Andersen, M.; Wernig, E.-M.; Hörmann, N. G.; Buller, N.; Reuter, K.; Kunze-Liebhäuser, J. Self-activation of copper electrodes during CO electro-oxidation in alkaline electrolyte. _Nat. Catal._ 2020, 1–7


# Efficient and Robust Jet Tagging at the LHC with Knowledge Distillation

Ryan Liu
University of California, Berkeley
Berkeley, CA 94720
Abhijith Gandrakota
Fermi National Accelerator Laboratory
Batavia, IL 60510
Jennifer Ngadiuba
Fermi National Accelerator Laboratory
Batavia, IL 60510
Maria Spiropulu
California Institute of Technology
Pasadena, CA 91125
Jean-Roch Vlimant
California Institute of Technology
Pasadena, CA 91125

###### Abstract

The challenging environment of real-time data processing systems at the Large
Hadron Collider (LHC) strictly limits the computational complexity of
algorithms that can be deployed. For deep learning models, this implies that
only models with low computational complexity that have weak inductive bias
are feasible. To address this issue, we utilize knowledge distillation to
leverage both the performance of large models and the reduced computational
complexity of small ones. In this paper, we present an implementation of
knowledge distillation, demonstrating an overall boost in the student models’
performance for the task of classifying jets at the LHC. Furthermore, by using
a teacher model with a strong inductive bias of Lorentz symmetry, we show that
we can induce the same inductive bias in the student model which leads to
better robustness against arbitrary Lorentz boost.

##  1 Introduction

In the past decades, deep learning has transformed the data analysis workflow
at high-energy physics experiments, achieving state-of-the-art performance on
many challenging tasks such as jet tagging, charged particle tracking,
particle flow algorithm, and neutrino event reconstruction [1, 2, 3, 4]. Such
breakthroughs mostly follow from novel architectures and large-scale datasets.
In particular, architectural innovation involves both general-purposed models
adopted from deep-learning research (e.g. convolutional neural networks,
transformers, graph neural networks [5, 1, 6]) and custom models that
incorporate strong inductive biases [7, 8]. However, these paradigms are not
applicable to each stage of the experimental data processing workflow,
especially to online applications such as the hardware trigger systems [9, 10]
of particle detectors at the LHC that are characterized by strict latency,
power, and resource constraints. For such applications, a generic method that
can improve deep learning models’ performance and robustness without scaling
the computation in the production environment becomes crucial.

Since knowledge distillation (KD) was proposed by Hinton in 2015 [11], it has
been widely used in natural language processing applications [12, 13]. KD is a
technique to transfer the learned knowledge from a large and cumbersome
“teacher” model to a more efficient “student” model. This is done by replacing
the ground truth (hard targets) with teacher-predicted probabilities (soft
targets). Remarkably, KD can significantly boost the performance of the
student model without any modification to the architecture itself. In this
paper, we describe an implementation of KD for jet tagging at the LHC and
demonstrate that powerful deep-learning models can be deployed to early-stage
event selections.

##  2 Related Work

###  2.1 Group Invariant Neural Networks

Inductive bias is a crucial part of designing a neural network. In the context
of jet physics, we have two important symmetries that shall be taken into
consideration: permutation invariance and Lorentz group invariance. Whereas
the constituents of a jet come in no particular order, the prediction should
not change upon the exchange of two particles in the input point cloud. At the
same time, since the source of the jet can have high momentum and cause the
jet to be boosted by a Lorentz transformation, the prediction should not
change by an arbitrary Lorentz transformation as well.

To incorporate permutation symmetry, it was proven that all permutation
invariant functions can be written as follows [14]:

$f(X)=\rho(\sum_{x\in X}\phi(x))$ (1)

We can parameterize $\rho$ and $\phi$ using multi-layer perceptrons (MLPs) to
make the network learnable. The resulting architecture is called a deep set
model. To improve the expressiveness of the group invariant neural network, we
can also consider attention-based models such as Transformers or Graph
Attention Networks [15, 16].

As for Lorentz symmetry, it is typically harder to incorporate its inductive
bias. Existing solutions include restricting the message-passing mechanism of
graph neural networks to use only Lorentz scalars such as $\langle x,y\rangle$
or $\|x-y\|^{2}$ [7] and manipulating representations of $SO^{+}(1,3)$ by
tensor-product and Clebsch–Gordan decomposition [8].

###  2.2 Knowledge Distillation

One of the common problems of the models described in the previous section is
that most of them do not scale well, especially for models that require pair-
wise computation such as Transformers and Graph Neural Networks. Therefore, we
propose to use KD to leverage the power of these models while not adding to
the inference-time cost. KD transfers learned knowledge by replacing the hard
target normally used to train the small and fast student model with soft
targets from the large and cumbersome teacher model. To be more precise, the
KD loss is defined as follows:

$L_{KD}(q;p,y)=(1-\lambda)\mathcal{H}(y,q)+\lambda
D_{KL}(\tilde{p}\|\tilde{q})$ (2)

where $q$ is the student’s output probabilities, $p$ is the teacher’s output
probabilities, and $y$ is the ground truth label. Probabilities with
$\tilde{\cdot}$ are the distributions softened by temperature $T$, i.e.,
$p(x)=\frac{e^{s(x)}}{\sum_{x^{\prime}}e^{s(x^{\prime})}}$ means
$\tilde{p}(x)=\frac{e^{s(x)}/T}{\sum_{x^{\prime}}e^{s(x^{\prime})/T}}$. The
first term $H(q,q)$ is the cross entropy loss with the ground truth as the
target, and the second term $D_{KL}(\tilde{p}\|\tilde{q})$ is the
Kullback–Leibler divergence with the softened teacher as the target. It was
pointed out that the KD loss has two major effects, namely, it can inject the
class relationship prior to the student and rescale the gradients of samples
according to the teacher’s relative confidence [17]. Moreover, recent studies
showed that KD is capable of transferring inductive bias [18]. That is, if the
teacher has a strong inductive bias, the student can learn from the KD loss to
generalize the same way as the teacher does.

##  3 Experiments

###  3.1 Top-quark jets Tagging at the LHC

In high energy physics, a jet refers to a collimated shower of particles that
result from the decay and hadronization of quarks $q$ and gluons $g$. The goal
of jet tagging is to identify which mother particle originated the jet such to
distinguish interesting and rare signal events from highly frequent background
ones. In this paper, we study the effect of KD on a common benchmark dataset
for jet tagging, the top tagging dataset [19]. The dataset contains signal
top-quark jets and background light quark and gluon jets simulated with
Pythia8 [20] and passed through a simulation of a typical LHC particle
detector obtained with Delphes [21]. There are 1.2M training events, 400K
validation events, and 400K testing events. Each jet is described by its
constituents’ four momenta, together with a flag indicating if it is a signal
event.

###  3.2 Model Architecture

In this paper, we demonstrate the efficacy of KD on two student models, the
deep set model [14] and an MLP model. The inputs are the constituents’ momenta
in cylindrical coordinates $(p_{T},\eta,\phi,m)$. Here we appended the
particle’s true mass $m$ to the input such that the model can distinguish
between different particle species. Note that the set of inputs we used here
are also available in a typical real-time data processing system thanks to the
particle flow algorithms at the Phase 2 CMS detector [22]. For the deep set
model, we used a 3-layer MLP with a hidden size of $128$ to parameterize both
$\rho$ and $\phi$, together with batch normalization [23] applied after leaky
relu activations [24]. The aggregation follows the design of the Energy Flow
Network such that per-particle embeddings are aggregated according to their
$p_{T}$ to enforce IR-safety [25]. For the MLP model, we sort particles
according to their $p_{T}$. Any jet with more than 128 particles is trimmed to
128 particles while jets with less than 128 particles are zero-padded. This
gives $4\times 128=512$ input features in total. We then feed them into an MLP
of 3 layers with 512 hidden features each, again equipped with leaky relu and
batch normalization. For more information, the code of this work is available
in this Github repository

###  3.3 Experiment Setup

To understand if KD can improve student performance and transfer inductive
biases, we picked the LorentzNet [7] as the teacher model. LorentzNet has a
very strong inductive bias of Lorentz-transformations invariance and is one of
the best-performing models on the top-tagging dataset. We designed two
experiments; (1) we evaluated the efficacy of KD for both the deep set and the
MLP models and compared their performance to the ones trained from scratch. We
experimented with temperatures of $T\in\\{1,3,5\\}$. (2) To show how KD can
transfer inductive bias, we trained the deep set model on an augmented dataset
that was boosted by $\beta$ uniformly sampled from $[0,\beta_{max}]$ along the
$x$-axis and evaluated the models on the boosted jets test set. We also chose
$\lambda=1$ for this experiment to prevent the model from learning from the
ground truth labels. We hypothesized that the ability for KD to transfer
inductive biases would be further improved when the student is exposed to some
corrupted samples that correspond to the teacher’s inductive bias. For each
model, we trained for 100 epochs with AdamW optimizer and a StepLR learning
rate scheduler.

###  3.4 Results

We report the accuracy, the area under curve (AUC), and the background
rejection ($\mathrm{R}_{X\%}$) of these models on the test set. The background
rejection is defined as $1/\mathrm{FPR}$ when the $\mathrm{TPR}$ is fixed to
$X\%$. Furthermore, with an eye to deployment to real-time system, we measured
the number of floating point operations (FLOPs) with the package fvcore [26],
which computes a model’s FLOPs with torch jit tracing. For the first
experiment, the results are reported in Table 1. Both the deep set and MLP
models showed performance gains from knowledge distillation. Remarkably, the
overall accuracy improvement for the MLP model is 1.5%, together with about a
factor of two improvement in background rejection. At the same time, the deep
set model showed about a 25% improvement in background rejection. The results
demonstrate the effectiveness of knowledge distillation in jet tagging.

Table 1: Comparison between models trained from scratch and knowledge
distillation.

#params FLOPs Accuracy AUC $\mathrm{Rej}_{30\%}$ $\mathrm{Rej}_{50\%}$ DeepSet
from scratch 68.2K 1.67M 0.930 0.9808 747 219 DeepSet KD $T=1$ 0.932 0.9818
926 241 DeepSet KD $T=3$ 0.932 0.9819 970 255 DeepSet KD $T=5$ 0.932 0.9819
970 248 MLP from scratch 527K 529K 0.904 0.9663 256 82 MLP KD $T=1$ 0.914
0.9726 375 119 MLP KD $T=3$ 0.918 0.9751 483 144 MLP KD $T=5$ 0.919 0.9750 503
146 LorentzNet (teacher) 224K 339M 0.942 0.9868 2195 498

(a) ROC curve of Deep sets

(b) ROC curve of MLPs

Figure 1: ROC curves of the two deep set and MLP models on the test dataset.
The knowledge distillation models have their curves above the models trained
from scratch.

For the second experiment, in all three configurations,
$\beta_{max}\in\\{0,0.1,0.3\\}$, we can see in Figure 1 that the robustness of
the deep set models significantly improved, presumably due to the inductive
biases transferred.

Finally, we observed that knowledge distillation can prevent models from
overfitting the data (see Figure 2(b)). This is especially useful when the
sample size is small. We hypothesize that this is due to a more complicated
learning objective that forces the model to generalize instead of memorize.

(a) Lorentz invariance test for the deep set model111The curve here for
LorentzNet is different from Figure 3 in [7]. This is because LorentzNet
included two auxiliary inputs that represent the beam line. In their paper,
the authors also boosted the beams. However, we believe that boosting the
beams encodes the information about the boosting so we decided not to follow
the convention..

(b) Training hard target loss

Figure 2: (a) we can see that the robustness with respect to Lorentz boosts
improved with KD in all three augmentation configurations. (b) In the first
experiment, the validation loss of MLP models with KD converges instead of
increasing as in the case without KD.

##  4 Conclusion

In this paper, we showed that knowledge distillation can improve jet tagging
performance and robustness. Without adding any additional operation, we
realized a $25\%$ improvement in background rejection. Furthermore, we
demonstrated that the teacher’s inductive bias can be transferred and help the
student generalize better. We hope that this work can serve as a base for
future deployment of deep learning models to real-time event selection systems
at the LHC, bringing the power of deep learning to the frontier of
experimental high-energy physics.

##  5 Broader Impact Statement

We expect that this work will stimulate the research and further discussions
on:

  1. 1.

Knowledge distillation between different architectures. Knowledge distillation
enables us to leverage the recent breakthroughs in large models to enhance the
small models we have. As shown in this work, knowledge distillation can work
for models of different architectures.

  2. 2.

Knowledge distillation as a method to improve robustness. In the Lorentz
invariance test, we have seen that knowledge distillation is capable of
improving student’s robustness. This may suggest that we can train a
cumbersome teacher model on an augmented dataset and transfer the robustness
to the student by knowledge distillation without having the student trained on
the augmented dataset, avoiding the accuracy and robustness trade-off.

##  6 Limitations

There exist a few limitations of this work. Firstly, we have only shown the
efficacy of knowledge distillation for a classification problem. Whether it
would work for other types of tasks in high-energy physics remains unclear.
Secondly, training the teacher model may be very costly. In this work, we used
a pre-trained model as the teacher. However, if one would like to train a
teacher model from scratch the cost of training may be significant.

##  7 Acknowledgement

This work is listed in Fermilab Technical Publications as FERMILAB-
PUB-23-748-CMS. AG and JN are supported by Fermi Research Alliance, LLC under
Contract No. DE-AC02-07CH11359 with the Department of Energy (DOE), Office of
Science, Office of High Energy Physics. JN and RL are also supported by the
U.S. Department of Energy (DOE), Office of Science, Office of High Energy
Physics “Designing efficient edge AI with physics phenomena” Project (DE-
FOA-0002705). JN is also supported by the DOE Office of Science, Office of
Advanced Scientific Computing Research under the “Real-time Data Reduction
Codesign at the Extreme Edge for Science” Project (DE-FOA-0002501). This work
was supported in part by the AI2050 program at Schmidt Futures (Grant
G-23-64934).

## References

  * [1] Huilin Qu, Congqiao Li, and Sitian Qian.  Particle transformer for jet tagging.  In International Conference on Machine Learning, pages 18281–18292. PMLR, 2022.
  * [2] Sylvain Caillou, Paolo Calafiura, Steven Andrew Farrell, Xiangyang Ju, Daniel Thomas Murnane, Charline Rougier, Jan Stark, and Alexis Vallier.  ATLAS ITk Track Reconstruction with a GNN-based pipeline.  Technical report, CERN, Geneva, 2022.
  * [3] Joosep Pata, Javier Duarte, Farouk Mokhtar, Eric Wulff, Jieun Yoo, Jean-Roch Vlimant, Maurizio Pierini, and Maria Girone.  Machine learning for particle flow reconstruction at CMS.  Journal of Physics: Conference Series, 2438(1):012100, feb 2023\.
  * [4] Alexander Shmakov, Alejandro Yankelevich, Jianming Bian, and Pierre Baldi.  Interpretable joint event-particle reconstruction for neutrino physics at nova with sparse cnns and transformers, 2023.
  * [5] CMS collaboration et al.  Identification of heavy, energetic, hadronically decaying particles using machine-learning techniques.  Journal of Instrumentation, 2020.
  * [6] Jonathan Shlomi, Peter Battaglia, and Jean-Roch Vlimant.  Graph neural networks in particle physics.  Machine Learning: Science and Technology, 2(2):021001, 2020.
  * [7] Shiqi Gong, Qi Meng, Jue Zhang, Huilin Qu, Congqiao Li, Sitian Qian, Weitao Du, Zhi-Ming Ma, and Tie-Yan Liu.  An efficient lorentz equivariant graph neural network for jet tagging.  Journal of High Energy Physics, 2022(7):1–22, 2022.
  * [8] Alexander Bogatskiy, Brandon Anderson, Jan Offermann, Marwah Roussi, David Miller, and Risi Kondor.  Lorentz group equivariant neural network for particle physics.  In International Conference on Machine Learning, pages 992–1002. PMLR, 2020.
  * [9] Technical Design Report for the Phase-II Upgrade of the ATLAS TDAQ System.  Technical report, CERN, Geneva, 2017.
  * [10] The Phase-2 Upgrade of the CMS Level-1 Trigger.  Technical report, CERN, Geneva, 2020.  Final version.
  * [11] Geoffrey Hinton, Oriol Vinyals, and Jeff Dean.  Distilling the knowledge in a neural network.  arXiv preprint arXiv:1503.02531, 2015.
  * [12] Xiaoqi Jiao, Yichun Yin, Lifeng Shang, Xin Jiang, Xiao Chen, Linlin Li, Fang Wang, and Qun Liu.  Tinybert: Distilling bert for natural language understanding.  arXiv preprint arXiv:1909.10351, 2019.
  * [13] Victor Sanh, Lysandre Debut, Julien Chaumond, and Thomas Wolf.  Distilbert, a distilled version of bert: smaller, faster, cheaper and lighter.  arXiv preprint arXiv:1910.01108, 2019.
  * [14] Manzil Zaheer, Satwik Kottur, Siamak Ravanbakhsh, Barnabas Poczos, Ruslan Salakhutdinov, and Alexander Smola.  Deep sets, 2018.
  * [15] Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N Gomez, Łukasz Kaiser, and Illia Polosukhin.  Attention is all you need.  Advances in neural information processing systems, 30, 2017.
  * [16] Petar Veličković, Guillem Cucurull, Arantxa Casanova, Adriana Romero, Pietro Lio, and Yoshua Bengio.  Graph attention networks.  arXiv preprint arXiv:1710.10903, 2017.
  * [17] Jiaxi Tang, Rakesh Shivanna, Zhe Zhao, Dong Lin, Anima Singh, Ed H. Chi, and Sagar Jain.  Understanding and improving knowledge distillation, 2021.
  * [18] Samira Abnar, Mostafa Dehghani, and Willem Zuidema.  Transferring inductive biases through knowledge distillation.  arXiv preprint arXiv:2006.00555, 2020.
  * [19] Gregor Kasieczka, Tilman Plehn, Jennifer Thompson, and Michael Russel.  Top quark tagging reference dataset, March 2019.
  * [20] Torbjörn Sjöstrand, Stefan Ask, Jesper R. Christiansen, Richard Corke, Nishita Desai, Philip Ilten, Stephen Mrenna, Stefan Prestel, Christine O. Rasmussen, and Peter Z. Skands.  An introduction to pythia 8.2.  Computer Physics Communications, 191:159–177, 2015.
  * [21] J De Favereau, Christophe Delaere, Pavel Demin, Andrea Giammanco, Vincent Lemaitre, Alexandre Mertens, and Michele Selvaggi.  Delphes 3: a modular framework for fast simulation of a generic collider experiment.  Journal of High Energy Physics, 2014(2):1–26, 2014.
  * [22] The Phase-2 Upgrade of the CMS L1 Trigger Interim Technical Design Report.  Technical report, CERN, Geneva, 2017.  This is the CMS Interim TDR devoted to the upgrade of the CMS L1 trigger in view of the HL-LHC running, as approved by the LHCC.
  * [23] Sergey Ioffe and Christian Szegedy.  Batch normalization: Accelerating deep network training by reducing internal covariate shift, 2015.
  * [24] Bing Xu, Naiyan Wang, Tianqi Chen, and Mu Li.  Empirical evaluation of rectified activations in convolutional network, 2015.
  * [25] Patrick T Komiske, Eric M Metodiev, and Jesse Thaler.  Energy flow networks: deep sets for particle jets.  Journal of High Energy Physics, 2019(1):1–46, 2019.
  * [26] Facebook Research.  fvcore.  https://github.com/facebookresearch/fvcore, 2019.