\section{#1}\setcounter{equation}{0}}
\renewcommand{\theequation}{\thesection.\arabic{equation}}
\newcommand{\be}{\begin{equation}}
\newcommand{\ee}{\end{equation}}
\newcommand{\bea}{\begin{eqnarray}}
\newcommand{\eea}{\end{eqnarray}}
\newcommand{\non}{\nonumber}
\newcommand{\ba}{\begin{array}}
\newcommand{\ea}{\end{array}}

\newcommand{\1}{{\underline{1}}}
\newcommand{\2}{{\underline{2}}}
\newcommand{\au}{\underline{a}}
\newcommand{\bu}{\underline{b}}
\newcommand{\cu}{\underline{c}}
\newcommand{\du}{\underline{d}}
\newcommand{\iu}{\underline{i}}
\newcommand{\ju}{\underline{j}}
\newcommand{\ku}{\underline{k}}
\newcommand{\lu}{\underline{l}}
\newcommand{\Iu}{\underline{I}}
\newcommand{\Ju}{\underline{J}}

\newcommand{\dsR}{{\mathbb R}}
\newcommand{\dsC}{{\mathbb C}}
\newcommand{\dsQ}{{\mathbb Q}}
\newcommand{\dsZ}{{\mathbb Z}}
\newcommand{\dsN}{{\mathbb N}}

\def\dt#1{{\buildrel {\hbox{\LARGE .}} \over {#1}}}

\newcommand{\bm}[1]{\mbox{\boldmath$#1$}}

\def\double #1{#1{\hbox{\kern-2pt $#1$}}}

\newcommand{\hm}{{\hat{m}}}
\newcommand{\hn}{{\hat{n}}}
\newcommand{\hp}{{\hat{p}}}
\newcommand{\hq}{{\hat{q}}}
\newcommand{\ha}{{\hat{a}}}
\newcommand{\hb}{{\hat{b}}}
\newcommand{\hc}{{\hat{c}}}
\newcommand{\hd}{{\hat{d}}}
\newcommand{\he}{{\hat{e}}}
\newcommand{\hM}{{\hat{M}}}
\newcommand{\hN}{{\hat{N}}}
\newcommand{\hA}{{\hat{A}}}
\newcommand{\hB}{{\hat{B}}}
\newcommand{\hC}{{\hat{C}}}
\newcommand{\hD}{{\hat{D}}}
\newcommand{\hE}{{\hat{E}}}
\newcommand{\hF}{{\hat{F}}}

\newcommand{\hal}{{\hat{\a}}}
\newcommand{\hbe}{{\hat{\b}}}
\newcommand{\hga}{{\hat{\g}}}
\newcommand{\hde}{{\hat{\d}}}
\newcommand{\hrh}{{\hat{\rho}}}
\newcommand{\hta}{{\hat{\tau}}}

\newcommand{\gd}{{\dot\g}}
\newcommand{\dd}{{\dot\d}}

\newcommand{\ts}{{\tilde{\s}}}
\newcommand{\tom}{{\tilde{\o}}}

\renewcommand{\Bar}{\overline}

\newcommand{\alu}{{\underline{\a}}}
\newcommand{\beu}{{\underline{\b}}}
\newcommand{\gau}{{\underline{\g}}}
\newcommand{\deu}{{\underline{\d}}}

\newcommand{\halu}{{\underline{\hal}}}
\newcommand{\hbeu}{{\underline{\hbe}}}
\newcommand{\hgau}{{\underline{\hga}}}
\newcommand{\hdeu}{{\underline{\hde}}}
\newcommand{\hrhu}{{\underline{\hrh}}}

\newcommand{\CD}{{\nabla}}

\newcommand{\sba}{{\bar{\s}}}

\newcommand{\te}{{\theta}}
\newcommand{\teb}{{\bar{\theta}}}

\newcommand{\sSp}{\mathsf{Sp}}
\newcommand{\sSU}{\mathsf{SU}}
\newcommand{\sSL}{\mathsf{SL}}
\newcommand{\sGL}{\mathsf{GL}}
\newcommand{\sSO}{\mathsf{SO}}
\newcommand{\sU}{\mathsf{U}}
\newcommand{\sS}{\mathsf{S}}
\newcommand{\sPSU}{\mathsf{PSU}}
\newcommand{\sPSL}{\mathsf{PSL}}
\newcommand{\sOSp}{\mathsf{OSp}}
\newcommand{\sSpin}{\mathsf{Spin}}
\newcommand{\sMat}{\mathsf{Mat}}

\newcommand{\bsubeq}{\begin{subequations}}
\newcommand{\esubeq}{\end{subequations}}
\newcommand{\opl}{{\oplus}}
\newcommand{\omn}{{\ominus}}

\newcommand{\ul}{\underline}
\newcommand{\eps}{{\ve}}
\newcommand{\bsigma}{\bar{\sigma}}
\newcommand{\dalpha}{{\dot{\alpha}}}
\newcommand{\dbeta}{{\dot{\beta}}}
\newcommand{\dgamma}{{\dot{\gamma}}}
\newcommand{\ddelta}{{\dot{\delta}}}
\newcommand{\dmu}{{\dot{\mu}}}
\newcommand{\dnu}{{\dot{\nu}}}
\newcommand{\dphi}{{\dot{\phi}}}
\newcommand{\btheta}{{\bar\theta}}
\newcommand{\N}{{\mathcal N}}
\newcommand{\eol}{\notag \\}
\newcommand{\rd}{\mathrm d}

\newcommand{\qW}{{\mathbf W}}
\newcommand{\qG}{{\mathbf G}}
\newcommand{\qH}{{\mathbf H}}
\newcommand{\qPsi}{{\mathbf \Psi}}
\newcommand{\qV}{{\mathbf V}}
\newcommand{\qOmega}{{\mathbf \Omega}}
\newcommand{\FF}{\nabla}
\newcommand{\bFF}{\bar \nabla}
\newcommand{\qomega}{{\bm\omega}}
\newcommand{\BCD}{{\bar\cD}}
\newcommand{\bphi}{{\bar\phi}}
\newcommand{\bj}{{\bar j}}
\newcommand{\gD}{{\mathbb D}}

\newcommand{\compint}{{\int \rd^4x\, e\, }}
\newcommand{\HC}{{\mathrm{c.c.}}}
\newcommand{\veps}{\varepsilon}
\newcommand{\de}{\nabla}

\newcommand{\halpha}{{\hat \alpha}}
\newcommand{\hbeta}{{\hat \beta}}
\newcommand{\hgamma}{{\hat \gamma}}
\newcommand{\hdelta}{{\hat \delta}}
\newcommand{\hmu}{{\hat \mu}}
\newcommand{\bnabla}{{\bm\nabla}}

\numberwithin{equation}{section}

\newcommand{\RM}{R(M)}
\newcommand{\RD}{R(\mathbb D)}
\newcommand{\RN}{R(N)}
\newcommand{\RJ}{R(\cJ)}
\newcommand{\RS}{R(S)}
\newcommand{\RK}{R(K)}
\newcommand{\RP}{R(P)}

\begin{document}

\begin{titlepage}
\begin{flushright}
Nikhef-2013-022\\
July, 2013\\
\end{flushright}
\vspace{5mm}

\begin{center}
{\Large \bf
Nonlocal action for the super-Weyl anomalies: \\
A new representation}
\\
\end{center}

\begin{center}

{\bf
Daniel Butter${}^{a}$ and Sergei M. Kuzenko${}^{b}$} \\
\vspace{5mm}

\footnotesize{
${}^{a}${\it Nikhef Theory Group \\
Science Park 105, 1098 XG Amsterdam, The Netherlands}}
~\\
\texttt{dbutter@nikhef.nl}\\
\vspace{2mm}

\footnotesize{
${}^{b}${\it School of Physics M013, The University of Western Australia\\
35 Stirling Highway, Crawley W.A. 6009, Australia}}
~\\
\vspace{2mm}

\end{center}

\begin{abstract}
\baselineskip=14pt
Using the recently discovered $\cN=1$ supersymmetric extension
of the conformal fourth-order scalar operator (introduced
originally by Fradkin and Tseytlin and also known as the
``Paneitz operator'' or ``Riegert operator''),
we derive a new representation for the nonlocal action generating
the super-Weyl anomalies.
\end{abstract}

\vfill

\vfill
\end{titlepage}

\newpage
\renewcommand{\thefootnote}{\arabic{footnote}}
\setcounter{footnote}{0}

\section{Introduction}

The concept of super-Weyl transformations \cite{HT, Siegel:SC} in the framework of
the Wess-Zumino formulation \cite{WZ,GWZ} for 4D $\cN=1$ supergravity
is a generalization of the ordinary Weyl transformations.
The property of a locally supersymmetric theory to be super-Weyl invariant
can naturally be recast in terms of the Ferrara-Zumino supercurrent multiplet \cite{FZ}
which contains the energy-momentum tensor and the spinor supersymmetry current along
with some other components.  It is pertinent for our subsequent discussion
to elaborate on this point in some more detail.

The Ferrara-Zumino  multiplet is realized  in terms of two constrained superfields:
a real vector  $T_{a} $ known as the supercurrent  and  a covariantly chiral
scalar $T$,  $\bar \cD_\ad T =0$, called the supertrace
(also known as the anomaly multiplet).
The supercurrent is the source of supergravity  \cite{OS,FZ78,Siegel77}.
Let $S [\c, H, \vf , \bar \vf] $ be the action of matter superfields $\c^i$ coupled
to background supergravity. The supergravity multiplet is fully described in terms
of the corresponding prepotentials $H^m$, $\vf$ and $\bar \vf$ \cite{Siegel78,SG},
of which the gravitational superfield \cite{OS} $H^m$ is real and the
conformal compensator \cite{Siegel78} $\vf$ is chiral.
The supercurrent and the supertrace originate as variational derivatives of the matter action
with respect to the supergravity prepotential. Specifically, it turns out that
\begin{subequations} \label{1.1}
\bea
T_a &=& \frac{\D } {\D H^a} S [\c, H, \vf , \bar \vf]~, \qquad \\
T &=& \frac{\d } {\d \s} S [\c, H, \re^\s \vf , \re^{\bar \s} \bar \vf]\Big|_{\s=\bar \s=0}~,
\qquad \bar \cD_\ad \s=0~,
\label{supertrace}
\eea
\end{subequations}
where $\s$ is an arbitrary covariantly chiral scalar
superfield.
Here $\D/ \D H^{a}$ denotes a covariantized variational derivative with respect to the gravitational superfield, see \cite{GGRS,BK} for pedagogical reviews.
If the matter fields are chosen to obey their equations of motion, $\d S/ \d \c^i =0$,
the condition that the matter action is locally supersymmetric is expressed
as the conservation equation
\bea
\bar \cD^\ad T_{\a\ad} +\frac{2}{3} \cD_\a T = 0~.
\eea

A super-Weyl transformation \cite{HT} is  a local rescaling of the chiral
compensator \cite{Siegel:SC}, $\vf \to \vf' = \re^\s \vf$, accompanied by a certain transformation
of the matter superfields of the form $\c^i \to \c'{}^i = \re^{-d_{(+)} \s - d_{(-)} \bar \s}\c^i$,
for some parameters $d_{(\pm)} $
(with the gravitational superfield $H^m$ being super-Weyl inert).
If the matter fields in \eqref{1.1} are chosen to satisfy
the equations of motion $\d S/ \d \c^i =0$, it follows from \eqref{supertrace}
that $T$ determines the super-Weyl variation of the action.
Without imposing the matter equations of motion,
the super-Weyl variation of the action is
\bea
\d_\s S =
\int \rd^4 x \,\rd^2 \q \,\cE \, \s T ~+ ~{\rm c.c.} ~+~\int \d_\s\c^i  \cdot \frac{\d S}{\d \c^i} ~,
\eea
where $\cE= \vf^3$ denotes the chiral integration measure
(see, e.g.,  \cite{GGRS,BK} for more details).
If the theory is super-Weyl invariant, then $T=0$ for the on-shell matter fields.
In a rigid supersymmetric limit, when $H^m \to 0$ and $\vf \to 1$,
such a theory becomes superconformal.

In the quantum theory, integrating out the matter fields leads to the
effective action $\G  [ H, \vf , \bar \vf]$. The quantum supertrace is
\bea
\langle T\rangle =  \frac{\d } {\d \s} \G [ H, \re^\s \vf , \re^{\bar \s} \bar \vf]\Big|_{\s=\bar \s=0}~.
\label{1.4}
\eea
Generically, the super-Weyl symmetry is anomalous at the quantum level.
This means that, if one starts from a classically super-Weyl invariant theory,
then in general $\langle T\rangle $ turns out to be a non-zero local functional
of the supergravity prepotentials, which cannot  be eliminated by adding local
counterterms to the effective action.
According to the cohomological analysis given in \cite{BPT},
the general  form of $\langle T\rangle $
in classically super-Weyl invariant theories  is
\bea
8\p^2 \langle T\rangle= c \,W^{\a\b\g} W_{\a\b\g} - a \,\cG
 +\frac{1}{16}h\, (\bar \cD^2 - 4R ) \cD^2 R ~,
\label{anomaly}
\eea
where $\cG$ denotes the following chiral scalar:
\bea
\cG =  W^{\a\b\g} W_{\a\b\g} -\frac{1}{4} (\bar \cD^2 - 4R ) (G^aG_a +2R \bar R)~.
\eea
It turns out that the functional  \cite{FZ78} (with $E$ the full superspace integration measure)
\bea
\int \rd^4 x \,\rd^2 \q \,\cE \, \cG = \int \rd^4 x \,\rd^2 \q \,\cE \,W^{\a\b\g} W_{\a\b\g}
+ \int \rd^4 x \,\rd^2 \q \,\rd^2 \bar \q \,E \, (G^aG_a +2R \bar R)
\eea
is a topological invariant (see  \cite{Buchbinder:1988tj,BK}
for the  proof),  which is related to the difference
of the Gauss-Bonnet and Pontryagin invariants.
The coefficients $a$ and $c$ in \eqref{anomaly}
are of nontrivial significance in general superconformal
field theories, see \cite{AFGJ,Osborn,SchT,KSch} and references therein.
On the other hand, the coefficient $h$ may be completely arbitrary, for its value may be freely changed
by adding to the effective action
a finite local counterterm proportional to
\bea
\int \rd^4 x \,\rd^2 \q \,\rd^2 \bar \q \,E \, R\bar R~.
\label{RbarR}
\eea

The explicit calculation of the anomalous supertrace \eqref{anomaly}
for  the chiral and vector multiplets was carried out in \cite{BK86}
(see \cite{BK} for a review).
This calculation made use of the super-$a_2$ (or, equivalently, super-$b_4$)
coefficients computed by McArthur \cite{McA2} (for the chiral multiplet)
using the superspace normal coordinates
\cite{McA} and independently in \cite{BK86} using the superfield Schwinger-DeWitt
technique. Given the anomalous supertrace  \eqref{anomaly}, it becomes possible to look for
a nonlocal effective action that generates the anomaly.
Such an action was derived in \cite{BK88} (see \cite{BK} for a review).
This was done with the aid of
a composite chiral scalar $\O$, introduced originally in \cite{GNSZ},
which obeys the massless equation
\bea
( \cD^2 -4\bar R ) \O=0~, \qquad \bar \cD_\ad \O =0
\eea
and possesses the super-Weyl transformation law\footnote{Using the super-Weyl transformation
of the antichiral torsion $\bar R$,
$\bar R \to \bar R' = -\frac{1}{4} \re^{-2\bar \s}
(\cD^2 -4 \bar R) \re^\s$, one can see that $\O$ defines a super-Weyl transformation
such that $\bar R' = 0$ and $R'=0$.}
\bea
\O \to \O' = \re^{-\s} \O~.
\eea
The explicit expression for $\O$ is
\bea
\O = 1+ \frac{1}{4 \Box_+} ( \bar \cD^2 -4R) \bar R,
\eea
where $\Box_+$ denotes the chiral d'Alembertian defined by
$\Box_+ \f = \frac{1}{16} ( \bar \cD^2 -4R) ( \cD^2 -4 \bar R) \f$,
for any covariantly chiral scalar $\f$.

Conceptually, the chiral superfield $\O$ is a supersymmetric analogue
of the composite scalar field $\o= 1 + \frac{1}{6} (\Box - \frac{1}{6} \cR)^{-1} \cR$,
with $\cR$ the scalar curvature,
introduced by Fradkin and Vilkovisky \cite{FV}.    The scalar field $\o$ was used
by Fradkin and Tseytlin \cite{FT84} to integrate the ordinary Weyl (or conformal) anomalies
(see \cite{DDI,Duff,DSch,Duff94} and references therein).
So the approach of \cite{BK88} may be thought of as a supersymmetric extension
of that given in \cite{FT84}. In the case of the Weyl anomaly,
there exists an alternative method \cite{Riegert}
to construct a nonlocal action generating the anomaly.\footnote{Nonlocal actions
for Weyl anomalies in  higher dimensions
were studied in \cite{KMM,Deser96,Deser2000}.}
It makes use of the conformal fourth-order scalar operator
\begin{align} \label{D_0}
\Delta_0 = \Box \Box +  \mathcal{D}^a \big(
	2 \mathcal{R}_{ab} \,\mathcal{D}^b
	- \tfrac{2}{3} \mathcal{R} \,\mathcal{D}_a
	\big)
\end{align}
discovered by Fradkin and Tseytlin\footnote{Fradkin and
Tseytlin \cite{FT1982} also constructed conformal
operators $\Delta_s$ for fields of different spin.}
in 1981 \cite{FT1982}
and re-discovered
by Paneitz in 1983 \cite{Paneitz}.\footnote{In the mathematics
literature, this operator is known as the Paneitz operator.
In the first preprint version of our paper, it was called
``the Paneitz-Riegert operator'' since the same operator had appeared
in \cite{Riegert}  and we were not aware of the earlier work by Fradkin and Tseytlin
\cite{FT1982}.}

It was mentioned in \cite{BK88} that the technique of \cite{Riegert} was not ``directly applicable
for the anomaly integration in curved superspace.'' The reason for that was the non-existence
of a super-Weyl covariant
chiral d'Alembertian\footnote{Given a chiral  d'Alembertian $\stackrel{\frown}{\Box}_+$ such that
$\stackrel{\frown}{\Box}_+ \J$ is covariantly chiral for any covariantly chiral scalar  $\J$,
the differential part of $\stackrel{\frown}{\Box}_+$ is uniquely fixed \cite{BK86},
$\stackrel{\frown}{\Box}_+ = \cD^a\cD_a +\frac{1}{4}R\cD^2 +\ri G^a \cD_a +\frac{1}{4}(\cD^\a R)\cD_\a +\cP$,  where the scalar superfield $\cP$ is covariantly chiral. The only free parameter in $\stackrel{\frown}{\Box}_+ $
is the chiral scalar $\cP$.}
$\stackrel{\frown}{\Box}_+$ such that the functional
\bea
\int \rd^4 x \,\rd^2 \q \,\rd^2 \bar \q \,E \, \bar \F \stackrel{\frown}{\Box}_+ \J
=\frac{1}{16} \int \rd^4 x \,\rd^2 \q \,\rd^2 \bar \q \,E \,  (\bar \cD^2 \bar \F ) \cD^2 \J + \dots
\label{1.11}
\eea
is super-Weyl invariant for any super-Weyl inert chiral scalars $\F$ and $\J$
(the ellipsis on the right of \eqref{1.11} denotes terms with two or fewer spinor derivatives).
The component counterpart of the operator in \eqref{1.11} is fourth-order in
vector derivatives.

A way out has only recently been found in
Ref. \cite{BdeWKL} which provided the required supersymmetric extension of
the Fradkin-Tseytlin $\Delta_0$ operator (using the novel
superspace formulation \cite{ButterN=1}  for $\cN=1$ conformal  supergravity).
Here we will use this operator to
derive a new representation for the nonlocal action generating
the super-Weyl anomaly.

 This paper is organized as follows. In section 2,  we describe the properties
 of the supersymmetric Fradkin-Tseytlin operator. In section 3,
 we derive the nonlocal effective action generating the super-Weyl anomaly.
Concluding comments are given in section 4.
 A brief summary of the Wess-Zumino superspace geometry,
following the notation and conventions of \cite{BK}, is given in the appendix.

\section{The supersymmetric Fradkin-Tseytlin operator}

In this section we give a detailed derivation of the $\cN=1$ super-Weyl covariant
operator which was introduced in \cite{BdeWKL} as the
supersymmetric extension of the Fradkin-Tseytlin  operator \eqref{D_0}.

As is well known, the algebra of covariant derivatives \eqref{algebra}
is invariant under the infinitesimal super-Weyl transformation \cite{HT}
associated with a chiral parameter $\s$, $\bar \cD_\ad \s=0$, and its complex conjugate $\bar \s$,
\begin{subequations}
\label{superweyl}
\bea
\d_\s \cD_\a &=& ( \hf \s - {\bar \s} )  \cD_\a - (\cD^\b \s) \, M_{\a \b}  ~, \\
\d_\s \bar \cD_\ad & = & ( \hf {\bar \s} - \s  )
\bar \cD_\ad -  ( \bar \cD^\bd  {\bar \s} )  {\bar M}_{\ad \bd} ~,\\
\d_\s \cD_{\a\ad} &=& -\hf( \s +\bar \s) \cD_{\a\ad}
-\frac{\ri}{2} (\bar \cD_\ad \bar \s) \cD_\a - \frac{\ri}{2} ( \cD_\a  \s) \bar \cD_\ad \non \\
&& - (\cD^\b{}_\ad \s) M_{\a\b} - (\cD_\a{}^\bd \bar \s) \bar M_{\ad \bd}~,
\eea
\end{subequations}
provided the torsion tensors transform\footnote{The super-Weyl transformation of $G_{\a\ad}$
given in \cite{BK}, eq. (5.5.14), contains a typo.}
 as follows:
\begin{subequations}
\label{superweyl-torsion}
\bea
\d_\s R &=& -2\s R -\frac{1}{4} (\bar \cD^2 -4R ) \bar \s
= (\bar \s - 2\s) R -\frac{1}{4} \bar \cD^2 \bar \s
~, \\
\d_\s G_{\a\ad} &=& -\hf (\s +\bar \s) G_{\a\ad} +\ri \,\cD_{\a\ad} (\bar \s- \s) ~,
\label{s-WeylG}\\
\d_\s W_{\a\b\g} &=&-\frac{3}{2} \s W_{\a\b\g}~.
\eea
\end{subequations}

Let $\bar \F$ be a covariantly antichiral scalar, $\cD_\a \bar \F=0$,
invariant under the super-Weyl transformations, $\d_\s \bar \F =0$.
It is an instructive exercise to check, using the algebra of covariant derivatives given in the Appendix,
the following super-Weyl  transformation law:
\bea
\d_\s \Big\{ \cD^2 \bar \cD^2 \bar \F + 8 \cD^\a (G_{\a\ad}\bar \cD^\ad \bar \F)\Big\}
&=& -(\s +\bar \s)  \Big\{ \cD^2 \bar \cD^2 \bar \F
+ 8 \cD^\a (G_{\a\ad} \bar \cD^\ad\bar \F)\Big\}
\non \\
&& -2 \bar \cD_\ad \Big\{ (\bar \cD^\ad \bar \F) \cD^2 \s
+ 4\ri (\cD^{\a\ad} \bar \F) \cD_\a \s \Big\} ~.
\label{second-order-operator}
\eea
To simplify the calculation, it is advantageous to make use of the identity
\bea
\d_\s \cD^2 = (\s-2\bar \s) \cD^2 + 2(\cD^\a \s) \cD_\a + \dots~,
\eea
where the ellipsis stands for all terms which involve the Lorentz generators
and annihilate  scalars.
It follows from \eqref{second-order-operator} that the operator $\D$ defined by
\bea
\D\bar \F := -\frac{1}{64} (\bar \cD^2 -4R ) \Big\{ \cD^2 \bar \cD^2 \bar \F
+ 8 \cD^\a (G_{\a\ad}\bar \cD^\ad \bar \F)\Big\}~, \qquad
\bar \cD_\ad \D\bar \F =0
\label{Delta}
\eea
has the super-Weyl transformation law
\bea
\d_\s \D\bar \F = -3\s \D\bar \F ~.
\label{sW-Delta}
\eea
We also introduce the conjugate operator
\bea
\bar\D \F := -\frac{1}{64} ( \cD^2 -4\bar R ) \Big\{ \bar \cD^2  \cD^2  \F
- 8 \bar \cD^\ad (G_{\a\ad} \cD^\a  \F)\Big\}~, \qquad
 \cD_\a \bar \D  \F =0~.
\eea
For any covariantly chiral scalars $\F$ and $\J$, it holds that
\bea
\int \rd^4 x \,\rd^2 \q \,\cE \, \J \D \bar \F &=&
\int \rd^4 x \,\rd^2 \bar \q\, \bar \cE \, \bar \F \bar \D \J \non \\
&=& \frac{1}{16}\int \rd^4 x \,\rd^2 \q \,\rd^2 \bar \q \,E \,
 \Big\{ (\cD^2 \J) \bar \cD^2 \bar \F
-8 (\cD^\a \J) G_{\a\ad}\bar \cD^\ad \bar \F\Big\} ~.~~~~
\label{symmetric}
\eea
In accordance with \eqref{sW-Delta},
this functional is super-Weyl invariant provided the chiral scalars $\F$ and $\J$
are inert under the super-Weyl transformations.
In the case that $\F = \J$, the functional \eqref{symmetric} is real.
Its component expression was given in \cite{FT} up to terms involving the gravitino.

The functional \eqref{symmetric} bears some resemblance to
\eqref{1.11} since both possess the same flat space limit, but the operators
underlying these two constructions are quite different. The chiral
d'Alembertian, which can be written
$\stackrel{\frown}{\Box}_+ = \frac{1}{16} (\bar \cD^2 - 4 R) \cD^2 + \cP$,
is a dimension-2 operator mapping a chiral multiplet to another chiral multiplet;
however, it cannot be made super-Weyl covariant for any choice of the chiral
function $\cP$, which was the point made in \cite{BK88}.
In contrast, the operator $\Delta$ underlying \eqref{symmetric}
is a dimension-3 super-Weyl covariant operator; it acts on
a weight-zero antichiral to yield a weight-three chiral superfield.
The overall chiral projector in \eqref{Delta} is critical for achieving
this \emph{manifest} super-Weyl covariance, because the quantity in
braces in \eqref{Delta} is neither super-Weyl covariant nor of definite
chirality by itself.
The resulting super-Weyl covariance of the functional \eqref{symmetric} is
absolutely critical for lifting the Fradkin-Tseytlin operator to superspace
and for enabling the construction which follows.

\section{Nonlocal effective action}

We now turn to deriving a nonlocal action that generates the super-Weyl anomaly \eqref{anomaly}.
We begin by introducing several important building blocks.

First, we require the super-Weyl transformations of some composite
expressions involving the torsion superfields. Using eq. \eqref{superweyl-torsion},
one can show
\bea
\d_\s (G^aG_a + 2R \bar R) &=& -(\s+\bar \s)  (G^aG_a + 2R \bar R)
-\hf \bar \cD_\ad (\bar R \bar \cD^\ad \bar \s)  - \hf \cD^\a (R\cD_\a \s) \non \\
&&-\hf \bar \cD^\ad (G_{\a\ad}   \cD^\a \s)+\hf  \cD^\a (G_{\a\ad}  \bar  \cD^\ad \bar \s)~,
\eea
which guarantees that the functional
\bea
\int \rd^4 x \,\rd^2 \q \,\rd^2 \bar \q\,  E \,(G^aG_a + 2R \bar R)
\eea
is super-Weyl invariant \cite{HT}.
Making use of the super-Weyl variation
\bea
\d_\s (\cD^2 R) = -(\s+\bar \s) \cD^2 R - 2\cD^\a(R\cD_\a \s) -\frac{1}{4} \cD^2\bar \cD^2 \bar \s~,
\eea
we further observe that
\bea
\d_\s \Big\{ G^aG_a + 2R \bar R  &-& \frac{1}{4} \cD^2 R    \Big\}  =
-(\s+\bar \s)  \Big\{ G^aG_a + 2R \bar R  - \frac{1}{4} \cD^2 R  \Big\}  \non \\
&+&\frac{1}{16} \cD^2\bar \cD^2 \bar \s
+\hf  \cD^\a (G_{\a\ad}  \bar  \cD^\ad \bar \s)
+\hf \bar \cD^\ad \Big( \bar R \bar \cD_\ad \bar \s
- G_{\a\ad}   \cD^\a \s\Big)~.
\eea
Now we introduce the composite chiral scalar
\bea
\X:= -\frac{1}{4} (\bar \cD^2 -4R ) \Big\{ G^aG_a + 2R \bar R
- \frac{1}{4} \cD^2 R \Big\}
\eea
and verify that
\bea
&& \d_\s \X = -3 \s \X + \D \bar \s~.
\eea

As a next step, we introduce two scalar Green's functions $G_{+-} (z, z') $ and $G_{-+} (z, z') $
that are related to each other by the rule
\bea
G_{+-} (z, z') = G_{-+} (z', z)
\eea
and obey the following conditions: \\
(i)   the two-point function $G_{-+} (z, z') $ is covariantly antichiral
in $z$ and chiral in $z'$,
\bea
\cD_\a  G_{-+} (z, z') =0~, \qquad \bar \cD_\ad' G_{-+} (z, z') =0~;
\eea
(ii) the two-point function $G_{-+} (z, z') $ satisfies the differential equation
\bea
\D G_{-+} (z, z') = \d_+ (z,z')~ .
\label{Green}
\eea
Here we have used the  chiral delta-function
\bea
 \d_+ (z,z'):= -\frac{1}{4} (\bar \cD^2 -4R ) E^{-1} \d^4(x-x')
 \d^2 (\q-\q') \d^2 (\bar \q -\bar \q') ~,
 \eea
which is covariantly chiral with respect to each of its arguments,
\bea
\bar \cD_\ad  \d_+ (z,z')=0~, \qquad \bar \cD'_\ad  \d_+ (z,z')=0~.
\eea
Under a finite super-Weyl transformation,  the delta-function
$ \d_+ (z,z')$ can be seen to change as
\bea
 \d_+ (z,z') \to \re^{-3\s}  \d_+ (z,z')~.
\label{sW-delta-function}
 \eea
It follows from the relations \eqref{sW-Delta}, \eqref{Green}
 and \eqref{sW-delta-function} that the Green's function $G_{-+} (z, z') $
is super-Weyl invariant, as is $G_{+-} (z, z') $.

Finally, it is useful to introduce a condensed notation
for the chiral and antichiral integration measures:
\bea
\int \rd \m_+ := \int \rd^4 x \,\rd^2 \q \,\cE ~, \qquad
\int \rd \m_- := \int \rd^4 x \,\rd^2 \bar \q\, \bar \cE ~.
\eea

With the building blocks given above, let us consider two nonlocal real functionals
\bea
\cF_1 &=&  \int \rd \m_+\int \rd \m'_-  \,\X(z) G_{+-} (z,z') \bar \X (z') ~,\\
\cF_2 &=&  \int \rd \m_+\int \rd \m'_- \, W^2(z) G_{+-} (z,z') \bar \X  (z')
~+ ~{\rm c.c.} ~,
\eea
with $W^2 = W^{\a\b\g}W_{\a\b\g} $.
The super-Weyl variations of these functionals are
\bea
\d_\s \cF_1 &=&  \int \rd \m_+ \, \s \X ~+ ~{\rm c.c.} ~, \\
\d_\s \cF_2 &=&  \int \rd \m_+ \, \s W^2 ~+ ~{\rm c.c.}
\eea
We also note that the super-Weyl variation of the local functional \eqref{RbarR} is
\bea
\d_\s \int \rd^4 x \,\rd^2 \q \,\rd^2 \bar \q \,E \, R \bar R =
\frac{1}{16} \int \rd \m_+ \, \s
(\bar \cD^2 -4R ) \cD^2 R ~+ ~{\rm c.c.}
\label{3.18}
\eea
It follows from these relations that the  effective action can be chosen in the form:
\bea
8\p^2 \G &=&  (c-a)  \int \rd \m_+\int \rd \m'_-  \, W^2(z)
G_{+-} (z,z') \bar \X  (z')   ~+ ~{\rm c.c.}\non \\
&-& a \int \rd \m_+\int \rd \m'_-  \, \X(z)
 G_{+-} (z,z') \bar \X  (z')
+ (h+a) \int \rd^4 x \,\rd^2 \q \,\rd^2 \bar \q \,E \, R \bar R ~.~~~
\label{3.19}
\eea

\section{Concluding comments}

The nonlocal action \eqref{3.19} is one of the main results of our paper.
It provides the generalization of Riegert's construction\footnote{Several papers
\cite{EO,Deser96,Deser2000} pointed out certain pathological properties
of the nonlocal action derived in \cite{Riegert}. More recently, this issue has been reconsidered in Ref.
\cite{Deser-rev} which ``rehabilitated'' \cite{Riegert}. Here we are only interested in the formal functional structure of the effective action.}
 \cite{Riegert} to $\cN=1$ supergravity.
The crucial property of $\G$ is that the corresponding supertrace $\langle T \rangle$,
eq. \eqref{1.4}, coincides with the super-Weyl anomaly \eqref{anomaly}.
The same anomaly is generated by the effective action derived twenty five
years ago in \cite{BK88} (see eq. (14) in \cite{BK88}).
The effective action derived in \cite{BK88}
differs from \eqref{3.19} by a super-Weyl invariant contribution
that depends only on the gravitational superfield and therefore vanishes in any conformally flat superspace.

In general, the super-Weyl anomaly may include not only the purely supergravity sector
\eqref{anomaly}, but also an additional contribution coming from a background vector supermultiplet, which is proportional to $\tr (W^\a W_\a)$, with $W_\a$ the covariantly chiral
field strength. Such a contribution is taken into account by an additional correction
to the effective action \eqref{3.19} obtained by the replacement
\bea
W^{\a\b\g} W_{\a\b\g}~ \to ~\tr (W^\a W_\a)~.
\eea

The functional \eqref{symmetric} is super-Weyl invariant for arbitrary
super-Weyl inert chiral scalars $\F$ and $\J$.  In a locally supersymmetric theory,
these chiral scalars may be chosen to be some of the dynamical variables or composite objects.
Thus \eqref{symmetric} generates a nontrivial higher derivative term in $\cN=1$ supergravity.

Because the functional \eqref{symmetric} is super-Weyl invariant, it is independent of the chiral compensator $\vf$ and its conjugate.  Therefore, this functional can naturally be
described within any off-shell version of $\cN=1$ supergravity.
As is known, any off-shell formulation for $\cN=1$ supergravity
can be realized by coupling a certain compensator to
$\rm U(1)$ superspace \cite{Howe} (see \cite{GGRS}
for a review), which gauges not only the Lorentz
group but the $\rm U(1)$ $R$-symmetry group as well.\footnote{In its turn, U(1) superspace \cite{Howe}
is a gauged-fixed version of  conformal superspace \cite{ButterN=1}.}
Using the same conventions
as in \cite{BK:Dual}, one finds that the operator $\Delta$ in $\rm U(1)$ superspace
is given by
\bea
\D\bar \F := -\frac{1}{64} (\bar \cD^2 -4R ) \Big\{ \cD^2 \bar \cD^2 \bar \F
+ 8 \cD^\a (G_{\a\ad}\bar \cD^\ad \bar \F)
+ \frac{16}{3} \bar X_\ad \bar\cD^\ad \bar \F\Big\}~,
\quad \bar \cD_\ad \D\bar \F =0~,~~~~~
\label{DeltaU(1)}
\eea
and obeys the super-Weyl transformation
\bea
\d_\L \D\bar \F = -3\L \D\bar \F ~,
\label{sW-DeltaU(1)}
\eea
where $\L$ is a real unconstrained super-Weyl transformation parameter.\footnote{We have flipped the
sign of $\L$ relative to \cite{BK:Dual}.}
In U(1) superspace the functional
\eqref{symmetric} turns into
\begin{align}
&\int \rd^4 x \,\rd^2 \q \,\cE \, \J \D \bar \F =
\int \rd^4 x \,\rd^2 \bar \q\, \bar \cE \, \bar \F \bar \D \J \non \\
& \quad =
\frac{1}{16}\int \rd^4 x \,\rd^2 \q \,\rd^2 \bar \q \,E \,
 \Big\{ (\cD^2 \J) \bar \cD^2 \bar \F
-8 (\cD^\a \J) G_{\a\ad}\bar \cD^\ad \bar \F
-\frac{16}{3} (\bar\cD_\ad \bar X^\ad) \J \bar\F \Big\} ~.~~~~
\label{symmetricU(1)}
\end{align}
This functional
is real when $\J = \F$, keeping in mind that
$\bar\cD_\ad \bar X^\ad = \cD^\a X_\a$. We note that the superspace
formulation of new minimal supergravity
can be derived from $\rm U(1)$ superspace by taking $R=0$,
so the expressions \eqref{DeltaU(1)} and \eqref{symmetricU(1)} hold
equally there.

As we already emphasized, the functional \eqref{symmetric} is super-Weyl invariant
under the condition that $\F$ and $\J$ are super-Weyl inert chiral scalars.
Choosing different types of superfields makes it possible to construct alternative
superconformal invariants. For instance, if $U$ and $V$ are super-Weyl inert
real scalars, then the functional
\bea
 \frac{1}{16}\int \rd^4 x \,\rd^2 \q \,\rd^2 \bar \q \,E \, U \cD^\a(\bar \cD^2 -4R)\cD_\a V
\eea
is super-Weyl invariant. In the case that $U=V$ is a dynamical superfield,
this functional coincides with the action of an Abelian vector multiplet \cite{WZ}.
On the other hand,  $U$ and $V$ may be composite, for instance
$U=V= K(\f^I, \bar \f^{\bar J}) $, where $K$ is the K\"ahler potential of a K\"ahler manifold,
and $\f^I$ are super-Weyl inert chiral scalars. Then, the above functional generates
a higher-derivative coupling which is invariant under arbitrary K\"ahler transformations
$K(\f, \bar \f) \to K(\f , \bar \f)  +\L (\f) +\bar \L(\bar \f)$, as well as under arbitrary
holomorphic isometries of the target K\"ahler  space.

The attempt to construct an $\cN=1$ supersymmetric extension
of the Fradkin-Tseytlin
operator was made a few years ago by Grosse \cite{Grosse}.
Unfortunately, the attempt was unsuccessful\footnote{One of us (SMK)
had incorrectly argued to Johannes Grosse in 2006 and Yu Nakayama in 2012
that a superfield version of $\Delta_0$ could not exist, based
on the non-existence of the super-Weyl covariant
d'Alembertian $\stackrel{\frown}{\Box}_+$. While working on \cite{BdeWKL},
it became apparent that the higher dimension operator $\Delta$ could
play the required role in deriving  a supersymmetric extension of  Riegert's construction.}
and this author concluded
that ``there is no superfield version of the Riegert operator for chiral fields of
Weyl weight 0.''
Although he did consider the integrand in the second line of
\eqref{symmetric}, he was not able to prove its super-Weyl invariance.
Most likely, this is due to its rather complicated transformation law.
We only discovered its invariance after first considering the associated
higher-order operator $\D$ defined by \eqref{Delta}, which proves to
possess the remarkably simple transformation law \eqref{sW-Delta}.
As demonstrated in \cite{BdeWKL}, the most natural origin of $\D$ is
within conformal superspace \cite{ButterN=1}.

In this paper we only considered the Wess-Zumino formulation for
$\cN=1$ supergravity \cite{WZ} which is characterized by  the old minimal set of auxiliary fields \cite{old}.
As is well known, there exist two other off-shell formulations for $\cN=1$ supergravity,
the non-minimal \cite{non-min,SG} and the new minimal \cite{new}.
In the case of non-minimal  supergravity, the super-Weyl transformations were
described by Siegel \cite{Siegel:SC}. For $n=-1$ non-minimal supergravity,
a superfield Fradkin-Tseytlin operator was constructed in \cite{Manvelyan}.
Super-Weyl transformations and trace anomalies in various off-shell versions
for $\cN=1$ supergravity have recently been discussed in \cite{Bonora}.
\\

\noindent
{\bf Acknowledgements:}\\
We thank Arkady Tseytlin for bringing the important references
\cite{FT1982} to our attention.
SMK is grateful to Stefan Theisen for a useful discussion and for pointing Ref. \cite{Deser-rev} out.
The work of DB was supported by
ERC Advanced Grant No. 246974, ``{\it Supersymmetry: a window to
non-perturbative physics}.''
The work of SMK   was supported in part by the Australian Research Council,
project No. DP1096372.


\section{Introduction}
At present,  QCD can not describe  quantitatively the low-energy
nucleon-antinucleon interaction, and various  phenomenological
approaches have been suggested in order to explain numerous
experimental data , see , e.g., Refs.
\cite{pot82,bonn91,bonn95,paris82,paris94,Partial94} and recent
reviews \cite{KBMR02,KBR05}. However, parameters of the models still
can not be extracted with a good accuracy  from the experimental
data \cite{Rich95}.

Very recently, renewed  interest in low-energy nucleon-antinucleon
physics has been stimulated by the experimental observation of a
strong enhancement of decay probability at low invariant mass of
$p\bar p$ in the processes  $J/\Psi\to \gamma p\bar p$ \cite{BES03},
$B^+\to K^+p\bar p$ and $B^0\to D^0 p\bar p$
\cite{Belle02a,Belle02b,BaBar05}, $B^+\to \pi^+ p\bar p$ and $B^+\to
K^0 p\bar p$ \cite{Belle04c}, $\Upsilon\to \gamma p\bar p$
\cite{CLEO06}. One of the most natural explanation of this
enhancement is final state interaction of the proton and antiproton
\cite{Ker04,Bugg04,Zou04,Loi05,Sib05,Sib06}.

A similar phenomenon was observed  in the investigation of the
proton (antiproton) electric, $G_E(Q^2)$, and  magnetic, $G_M(Q^2)$,
form factors  in the process $e^+e^- \to p\bar p$
\cite{Bardin94,Armstrong93,Aubert06}. Namely, it was found that the
ratio $|G_E(Q^2)/G_M(Q^2)|$ strongly depends on $Q^2=4E^2$ (in the
center-of-mass frame) in the narrow  region of  the energy $E$  near
the  threshold of $p\bar p$ production. Such  strong dependence at
small $E$ is related to the large-scale interaction of proton and
antiproton. Therefore, it is possible to apply the approaches of
\cite{pot82,bonn91,bonn95,paris82,paris94,Partial94} for an
explanation of experimental data in the process  $e^+e^- \rightarrow
p\bar p$. In the present paper, we use the Paris nucleon-antinucleon
optical potential $V_{N\bar N}$  which has the form
\cite{paris82,paris94}:
\begin{equation}\label{VNN}
V_{N\bar N} = U_{N\bar N}-i\, W_{N\bar N} \, ,
\end{equation}
where the real part $U_{N\bar N}$  is the  $G$-parity transform of
the well established Paris $NN$ potential for the long- and
medium-range distances ($r\gtrsim 1\mbox{fm}$), and some
phenomenological part for the short distances. The absorptive part
$W_{N\bar N}$ of the optical potential takes into account the
inelastic channels of $N\bar N$ interaction, i.e. annihilation into
mesons. It is essential at short distances and depends on the
kinetic energy of the particles. We perform calculations in the non-relativistic approximation.
The Coulomb interaction between proton and antiproton is important only for the kinetic energy  $T\lesssim (\pi\alpha)^2 M \sim 1$~MeV, where $\alpha$ is the fine structure constant, and $M$ is the proton mass. Here we consider the process for kinetic energies $T\gg 1$~MeV, therefore we neglect the Coulomb interaction.

Taking into account that the difference of the cross sections  $e^+e^- \rightarrow
p\bar p$ and  $e^+e^- \rightarrow n\bar n$ is small \cite{Aubert06,Antonelli98},  we calculate  the cross sections at a given isospin final states and   compare them  with the
 experimental data for $e^+e^- \rightarrow p\bar p$. As a result we found that the amplitude with the isospin $I=0$  strongly dominates.
  Besides, our prediction for the ratio  $|G_E(Q^2)/G_M(Q^2)|$ depends on the parameters of the Paris potential but independent
  of the form factor at threshold. The ratio is also in a qualitative agreement with experimental data.

\section{Amplitude of the process}
In the nonrelativistic  approximation, the amplitude of $N\bar{N}$
pair production in a certain isospin channel $I=0,\,1$ near  threshold can be presented as follows (in units $4\pi\alpha/Q^2$):
\begin{eqnarray}\label{2}
&&T_{\lambda\mu}^I=\sqrt{2} \bm{\epsilon}_\lambda^* \Bigg[\mathcal
G_s^I\bm e_\mu +\mathcal G_d^I \frac{{\bm k}^2{\bm e}_\mu-3({\bm k}\cdot{\bm e}_\mu){\bm k}}{6M^2}\Bigg]\, ,\nonumber\\
&&\mathcal G_s^I=\mathcal F_1^I(Q^2)+\mathcal F_2^I(Q^2)+\frac{\beta^2}{6}[\mathcal F_2^I(Q^2)-\mathcal F_1^I(Q^2)]\,,\nonumber\\
&&\mathcal G_d^I=\mathcal F_1^I(Q^2)-\mathcal F_2^I(Q^2)\,,
\end{eqnarray}
where $\beta=k/M\ll 1$ , ${\bf e}_\mu$ is a virtual photon polarization vector, corresponding to the projection
of spin  $J_z=\mu=\pm 1$,  and $\bm{\epsilon}_\lambda$ is the spin-1 function of $N\bar{N}$ pair, $\lambda=\pm1,\, 0$ is the projection of spin on the vector $\bm k$. Two tensor structures in Eq.(\ref{2}) correspond to the s-wave and
d-wave production amplitudes. The total angular momentum of the
$N\bar{N}$ pair is fixed by a production mechanism. The functions
$\mathcal F_1^I(Q^2)$ and $\mathcal F_2^I(Q^2)$ are the Dirac form factors of the $N\bar{N}$ pair which include the  effects of final state interaction and have  a pronounced $Q^2$ behavior near the threshold. Summation over the polarization of  nucleon  pair and averaging over the  polarization of   virtual photon is performed using the equations,
\begin{equation}
\sum_{\lambda=1,2,3} \epsilon_\lambda^{i*}\epsilon_\lambda^j=\delta^{ij}\,,\quad
 \frac{1}{2}\sum_{\mu=1,2} e_\mu^{i*}e_\mu^j=\frac{1}{2}\delta^{ij}_\perp=\frac{1}{2}(\delta^{ij}-P^iP^j/P^2)\,,
\end{equation}
where $\bm P$ is the electron  momentum.

Our aim is to single out the  effects of final state interaction. In order to do that, we write the amplitude (\ref{2}) in the form
\begin{eqnarray}\label{3}
T_{\lambda\mu}^I & =&\sqrt{2} \int \frac{d^3p}{(2\pi)^3}{\bm\Phi}^{I(-)*}_{\bm k\lambda}({\bm p})\cdot \Bigg[G_s^I\bm e_\mu +G_d^I\frac{{\bm p}^2{\bm e}_\mu-3({\bm p}\cdot{\bf
e}_\mu){\bm p}}{6M^2}\Bigg]\, ,
\end{eqnarray}
where $\bm\Phi^{I(-)}_{\bm k\lambda}({\bm p})$ is the Fourier transform of the function $\bm\Psi^{I(-)}_{\bm k\lambda}({\bm r})$, the wave function of the $N\bar{N}$ pair in  coordinate space.
This wave function is the  solution of the Schr\"odinger equation
\begin{equation}
\mathbf{\Psi}^{I(-)*}_{{\bm k}\lambda}({\bm r})\hat{H}=M\beta^2
\mathbf{\Psi}^{I(-)*}_{{\bm k}\lambda}({\bm r})\,,\quad
\hat{H}=\frac{{\bm p}^2}{M}+V_{N\bar{N}}\,,
\end{equation}
where $V_{N\bar{N}}$ is the optical potential. Note that
$\mathbf{\Psi}^{I(-)*}_{{\bm k}\lambda}({\bm r})$ is the left
eigenfunction of the bi-orthogonal set of eigenfunctions of the
non-Hermitian operator $\hat{H}$. The asymptotic form of the wave function at large distances reads
\begin{equation}
\bm\Psi^{I(-)}_{\bm k\lambda}({\bm r})\approx \bm{\epsilon}_\lambda\,e^{i\bm k\cdot\bm r}+f_{\lambda\lambda'} \frac{e^{-ikr}}{r}\bm{\epsilon}_{\lambda'} \,.
\end{equation}
In Eq.(\ref{3}) the form factors $G_s$ and $G_d$  are
\begin{eqnarray}\label{31}
&&G_s^I=F_1^I+ F_2^I+\frac{\beta^2}{6}(F_2^I-F_1^I)\,,\nonumber\\
&& G_d^I=F_1^I- F_2^I\,,
\end{eqnarray}
where the "bare" Dirac form factors  $F_1$ and  $F_2$  do not account for the effect of
final state interaction. Near threshold  these form factors  are  smooth
functions of $Q^2$ and can be treated  as  phenomenological constants.

\section{Wave function}

Let us introduce the vector spherical functions $\bm Y_{J\,\mu}^L(\bm n)$ as
\begin{eqnarray}\label{A}
\bm Y_{J\,\mu}^L(\bm n)=\sum_m  C^{J\,\mu}_{L\,m,\,1\,\mu-m}\,Y_{L\,m}(\bm n) \bm \epsilon_{\mu-m}\,,
\end{eqnarray}
where $Y_{L\,m}(\bm n)$ are spherical harmonics,  $C^{J\mu}_{L\,m,\,1\,\mu-m}$ are
Clebsch-Gordan coefficients, and $\bm n=\bm r/r$. In Eq. (\ref{A}) the quantization axes is directed along the vector $\bm k$. Then the wave function $\bm\Psi^{I(-)}_{\bm k\lambda}({\bm r})$ can be written in the form
\begin{eqnarray}\label{WFpsi}
&&\bm\Psi^{I(-)}_{\bm k\lambda}({\bm r})=\sum_J \sqrt{4\pi(2J+1)} C^{J\,\lambda}_{J\,0,\,1\,\lambda} v^{I*}_J(r)\bm Y_{J\,\lambda}^J(\bm n)\nonumber\\
&&+\sum_J \sqrt{4\pi(2J-1)} C^{J\,\lambda}_{J-1\,0,\,1\,\lambda}\Big[ u^{I*}_{1J}(r)\bm Y_{J\,\lambda}^{J-1}(\bm n)+w^{I*}_{1J}(r)\bm Y_{J\,\lambda}^{J+1}(\bm n)\Big]                              \nonumber\\
&&+\sum_J \sqrt{4\pi(2J+3)} C^{J\,\lambda}_{J+1\,0,\,1\,\lambda}\Big[ u^{I*}_{2J}(r)\bm Y_{J\,\lambda}^{J-1}(\bm n)+w^{I*}_{2J}(r)\bm Y_{J\,\lambda}^{J+1}(\bm n)\Big]\,,                               \end{eqnarray}
 Here the functions $v^I_J(r)$,  $u^I_{nJ}(r)$, and $w^I_{nJ}(r)$ have the asymptotic form at large distances
\begin{eqnarray}\label{asvuw}
&&v^I_J(r)=\frac{1}{2ikr}\Big[S_0^{IJ}e^{i[kr-J\pi/2]}-e^{-i[kr-J\pi/2]}\Big]\,,\nonumber\\
&&u^I_{1J}(r)=\frac{1}{2ikr}\Big[S_{11}^{IJ}e^{i[kr-(J-1)\pi/2]}-e^{-i[kr-(J-1)\pi/2]}\Big]\,,\nonumber\\
&&w^I_{1J}(r)=\frac{1}{2ikr}S_{12}^{IJ}e^{i[kr-(J+1)\pi/2]}\,,\nonumber\\
&&u^I_{2J}(r)=\frac{1}{2ikr}S_{21}^{IJ}e^{i[kr-(J-1)\pi/2]}\,,\nonumber\\
&&w^I_{2J}(r)=\frac{1}{2ikr}\Big[S_{22}^{IJ}e^{i[kr-(J+1)\pi/2]}-e^{-i[kr-(J+1)\pi/2]}\Big]\,,
\end{eqnarray}
where $S_0^{IJ}$  and $S_{ij}^{IJ}$ are some functions of energy with $|S_0^{IJ}|\leqslant 1$ and $|S_{ij}^{IJ}|\leqslant 1$. Due to angular momentum conservation, only the terms with $J=1$ and $L=J\pm 1$ (i.e., $L=0,\,2$) contribute to the matrix element (\ref{3}). Then the amplitude (\ref{3}) can be written as
\begin{eqnarray}\label{dress1}
T_{\lambda\mu}^I & =&\sqrt{2}\lim_{r\rightarrow 0}\Bigg[G_s^I\bm e_\mu  -G_d^I\frac{{\bm e}_\mu{\triangle}-3({\bm \nabla}\cdot{\bf e}_\mu){\bm \nabla}}{6M^2}\Bigg]\mathbf{\psi}^{I*}_{{\bm k}\lambda}(\bm r)\, ,\nonumber\\
&&\psi^{I}_{{\bm k}\lambda}(\bm r)= [u^{I*}_{11}(r)\bm \epsilon_\lambda+w^{I*}_{11}(r)\sqrt{4\pi}\bm Y_{1\,\lambda}^2(\bm n)]\nonumber\\
&&+ \sqrt{5}\, C^{1\,\lambda}_{2\,0,\,1\,\lambda}[u^{I*}_{21}(r)\bm \epsilon_\lambda+w^{I*}_{21}(r)\sqrt{4\pi}\bm Y_{1\,\lambda}^2(\bm n)]\,.
\end{eqnarray}
Finally we have
\begin{eqnarray}\label{21}
&&\mathcal G_s^I= G_s^I\,u^{I}_{11}(0)+\frac{5\, G_d^I }{\sqrt{2}M^2}\lim_{r\rightarrow 0}\left(\frac{w^{I}_{11}(r)}{r^2}\right)\,,\nonumber\\
&&\mathcal G_d^I=\frac{6G_s^I}{\sqrt{2}\beta^2}\,u^{I}_{21}(0)+15\, G_d^I \lim_{r\rightarrow 0}\left(\frac{w^{I}_{21}(r)}{k^2r^2}\right)\,.
\end{eqnarray}
The first term in $\mathcal G_d^I$ contains the large factor $6/\beta^2$ while the second term in $\mathcal G_s^I$ is small due to the proton mass $M$ in denominator. Thus, in the non-relativistic approximation the amplitude $T_{\lambda\mu}^I$ reads,
\begin{eqnarray}\label{dress2}
T_{\lambda\mu}^I & =&G_s^I\Bigg\{ \sqrt{2}u^{I}_{11}(0)(\bm e_\mu\cdot \bm{\epsilon}_\lambda^*)+
 u^{I}_{21}(0)[ ({\bm e}_\mu\cdot \bm{\epsilon}_\lambda^*)-3(\hat{\bm k}\cdot  {\bm e}_\mu)(\hat {\bm k}\cdot\bm{\epsilon}_\lambda^*)]\Bigg\}\, \,,
\end{eqnarray}
where $\hat{\bm k}=\bm k/k$. The interpretation of this equation is the following.  As a result of re-scattering due to the tensor forces, the pair  produced  at the origin in $s$-wave has a non-zero amplitude to transfer to $d$-wave.

\section{Cross section and Sachs form factors}

The cross section corresponding to the amplitude (\ref{3}) has the form in the center-of-mass frame (see, e.g., Ref.~\cite{BLP1971})
\begin{eqnarray}\label{CS0}
\frac{d\sigma}{d\Omega} & =&\frac{\beta\alpha^2}{4Q^2}\,\Bigg[|G_M(Q^2)|^2(1+\cos^2\theta)+\frac{4M^2}{Q^2}\,
|G_E(Q^2)|^2\sin^2\theta\Bigg]\,.
\end{eqnarray}
Here $\theta$ is the angle between the electron (positron) momentum $\bm P$ and the momentum of the final particle $\bm k$. In terms of the "dressed" form factors ${\cal G}_s^I$ and ${\cal G}_d^I$ the electromagnetic Sachs form factors, corresponding to the contribution of the amplitude with the isospin $I$,  have the form
\begin{eqnarray}\label{FFnr}
&&G_M^I={\cal G}_s^I+\frac{\beta^2}{6}{\cal G}_d^I=G_s^I[u^{I}_{11}(0)+\frac{1}{\sqrt{2}}u^{I}_{21}(0)]
\,,\nonumber\\
&& \frac{2M}{Q}G_E^I={\cal G}_s^I-\frac{\beta^2}{3}{\cal G}_d^I=G_s^I[u^{I}_{11}(0)-{\sqrt{2}}u^{I}_{21}(0)] \,.
\end{eqnarray}
Thus, in the non-relativistic approximation the ratio $G_E^I/G_M^I$ is independent of the constant $G_s^I$,
\begin{eqnarray}\label{GEGMnr}
&&\frac{G_E^I}{G_M^I}=
\dfrac{u^{I}_{11}(0)-{\sqrt{2}}u^{I}_{21}(0)}{u^{I}_{11}(0)+\dfrac{1}{\sqrt{2}}u^{I}_{21}(0)}
\,.
\end{eqnarray}
Note that the electromagnetic interaction is important only in the narrow region where $\beta\sim \pi\alpha$ and the nucleon energy is $E=M\beta^2/2\sim 0.3$MeV. In this paper we will not consider this narrow region and neglect the electromagnetic interaction in the potential. Then,
the amplitude of $p\bar{p}$ pair production, $T^{(p)}_{\lambda\mu}$,  and the amplitude of $n\bar{n}$ pair production,  $T^{(n)}_{\lambda\mu}$, has the form
$$T^{(p)}_{\lambda\mu}=\frac{T^{1}_{\lambda\mu}+T^{0}_{\lambda\mu}}{\sqrt{2}}\,,\quad
T^{(n)}_{\lambda\mu}=\frac{T^{1}_{\lambda\mu}-T^{0}_{\lambda\mu}}{\sqrt{2}}\,.$$
The contribution of the isospin $I$ to the total cross section of the nucleon pair production reads
\begin{eqnarray}\label{CStot}
\sigma^I & =&\frac{2\pi\beta\alpha^2}{Q^2}\,|G_s^I|^2\,[|u^I_{11}(0)|^2+|u^I_{21}(0)|^2]\,.
\end{eqnarray}
Thus, to describe the energy dependence of the ratio ${G_E^I}/{G_M^I}$ and the cross section
$\sigma^I$ in the non-relativistic approximation, it is necessary to know the functions
$u^I_{11}(0)$ and $u^I_{21}(0)$.

Let us write the hamiltonian  $H^I$ for the isospin $I$ as follows,

\begin{eqnarray}
&&H^I=\frac{p_r^2}{M}+V_0^I(r)\delta_{L0}+V_2^I(r)\delta_{L2}+V_3(r)\,S_{12}\,,\nonumber\\
&&S_{12}=6(\bm S\cdot\bm n)^2-4\,,
\end{eqnarray}
where $\bm S$ is the spin operator for the  spin-one system of produced pair,  $(-p_r^2)$ is the radial part of the Laplace operator, and $L$ denotes the orbital angular momentum.  Then the radial wave functions $u^I_{n1}$ and $w^I_{n1}$ , $n=1,\, 2$~, satisfy  the equations

\begin{eqnarray}\label{equation}
&&\frac{p_r^2}{M}\chi +{\cal V}\chi=2E\chi\,,\nonumber\\
&&{\cal V}=\begin{pmatrix}V_0^I & -2\sqrt{2}\,V_3^I \\
-2\sqrt{2}\,V_3^I  & V_2^I-2V_3^I
\end{pmatrix}\,,\quad \chi=\begin{pmatrix}u^{I}_{n1}\\w^{I}_{n1}\end{pmatrix}\,.
\end{eqnarray}
The asymptotic form of the solutions at large distances is, Eq. (\ref{asvuw}),
\begin{eqnarray}\label{asvuw1}
&&u^{I}_{11}(r)=\frac{1}{2ikr}\Big[S_{11}^{I1}\,e^{ikr}-e^{-ikr}\Big]\,,\nonumber\\
&&w^{I}_{11}(r)=-\frac{1}{2ikr}S_{12}^{I1}\,e^{ikr}\,,\nonumber\\
&&u^{I}_{21}(r)=\frac{1}{2ikr}S_{21}^{I1}\,e^{ikr}\,,\nonumber\\
&&w^{I}_{21}(r)=\frac{1}{2ikr}\Big[-S_{22}^{I1}e^{ikr}+e^{-ikr}\Big]\,.
\end{eqnarray}

\section{Results and discussion}

It is known that the difference of the cross section of the processes  $e^+e^- \rightarrow
p\bar p$ and  $e^+e^- \rightarrow n\bar n$ is small \cite{Aubert06,Antonelli98}. Therefore it is natural to suggest that one of the isospin amplitudes is much larger than  another one.
Fig.~1 shows the cross section $\sigma^0$ at $| G_s^0|^2=101.8$ (solid line) as well as the cross section  $e^+e^- \rightarrow p\bar p$. The only free parameter $| G_s^0|^2$ has been found by normalizing the theoretical curve to the data at the third experimental point above threshold.
It is seen that the cross section for $I=0$ perfectly reproduces the shape of the data in a wide range of energy near threshold (from threshold up to $2E=2.2$~GeV). Besides, for this case
the value of $| G_s^0|=|F_1^0+F_2^0|\approx 10$ looks reasonable compared to the value $ |G_s^1|\approx 100$ obtained by assumption of isovector dominance.

Fig.~2 shows the cross section $\sigma^1$  normalized at the same data point. It is seen that the shape of the theoretical curve does not reproduce the energy behaviour of the measured cross section. In addition, the value of the fitting parameter $ |G_s^1|^2=10184$ does not look reasonable.

Fig.~3 shows the ratio $|{G_E^0}/{G_M^0}|$, see Eq.(\ref{GEGMnr}), together with the experimental data Ref.\cite{Aubert06}. Our prediction is in a qualitative agreement with the data.

Isoscalar dominance in a model of final-state interaction based on the Paris potential has a simple explanation. At very small distances, the Paris potential corresponding to $I=1$ is strongly repulsive while the potential for  $I=0$ is strongly  attractive. As a result, the wave function for $I=1$ at small distances is strongly suppressed as compare to  the case   $I=0$.

To conclude, using the Paris optical potential, we show that the isoscalar amplitude dominates in the cross section of the process   $e^+e^- \rightarrow N\bar N$. Our prediction for the cross section is in good agreement with the data. The prediction for  $|{G_E^0}/{G_M^0}|$ obtained under assumption of the isoscalar dominance agrees qualitatively with the data. To confirm our statement on the isoscalar dominance, other optical potential models should be considered.

\section*{Acknowledgements}
The work  was supported by the Ministry of Education and Science of the
Russian Federation.


\section{Introduction}\label{sec.intro}

The goal of this paper is to introduce a class of new integrable models
with the help of exact results in supersymmetric gauge theories.
The integrability here refers to solutions of the Yang-Baxter equation
(YBE).

Our integrable model will be constructed from the following relation between
4d $\ensuremath{\mathcal{N}}=1$ quiver gauge theories and the 2d integrable models,
which we call the {\it Gauge/YBE correspondence}:
\begin{equation}\begin{aligned}
\ensuremath{\mathcal{I}}_\textrm{4d index}[S^1\times S^3/\ensuremath{\mathbb{Z}}_r]= Z_\textrm{2d spin}[C] \ .
\label{GaugeYBE}
\end{aligned}\end{equation}

On the left hand side of \eqref{GaugeYBE}
we have the 4d lens index \cite{Benini:2011nc}
for a class of 4d $\ensuremath{\mathcal{N}}=1$ supersymmetric $SU(N)$ quiver gauge theories
defined from quiver diagrams on a
two-dimensional geometry $C$, which is either
a disc \cite{Xie:2012mr,Franco:2012mm}
or a torus \cite{Hanany:1997tb,Hanany:1998it,Hanany:2005ve}.
The 4d lens index
is a twisted partition function on $S^1\times S^3/\ensuremath{\mathbb{Z}}_r$
and is a function of several fugacities.
Two of the fugacities will be denoted by $p,q$,
and others for global symmetries will be parametrized by
$R_i$.

One the right hand side of \eqref{GaugeYBE}
we have the
partition function of a 2d classical integrable spin system
defined on $C$, where the quiver diagram is identified with the lattice of the spin chain.
The boundary condition of the spin chain is either periodic or fixed
depending on whether $C$ is a torus or a disc.
The parameters of the index, $p,q$ and $R_i$, are translated into the
temperature-like parameters and the
spectral parameters of the spin chain, respectively.

The correspondence \eqref{GaugeYBE}
holds for each value of $r$ and $N$.
The $r=1$ case of \eqref{GaugeYBE} was discussed in the
earlier paper by the author
and the collaborators \cite{Yamazaki:2012cp,Terashima:2012cx,Xie:2012mr}.
There the 4d $S^1\times S^3$ index \cite{Kinney:2005ej,Romelsberger:2005eg}
of the 4d gauge theories
was identified with the
statistical mechanical integrable model discussed
in  \cite{Bazhanov:2010kz,Bazhanov:2011mz,Bazhanov:2013bh}
(see also \cite{Spiridonov:2010em}).

For the more general case $r>1$ discussed in this paper,
we can use the correspondence \eqref{GaugeYBE} to {\it define}
a new statistical model from the known gauge answer on the left hand side,
for each pair of positive integers $(r,N)$.
The spin at a vertex of the lattice has both discrete and
continuous spins, each of which has $N$ components with $1$ constraint imposed.
The R-matrix is written in terms of lens elliptic gamma functions
(Appendix \ref{sec.gamma}),
and is a function of
two temperature-like parameters $p,q$ as well as the spectral parameters
$R_i$.
Our solution to YBE is one of the most general known in the literature,
see more comments in section \ref{sec.comments}.

The origin of the integrability of the 2d spin system
has a clear-cut explanation
on the
supersymmetric gauge theory side of the
correspondence \eqref{GaugeYBE}:
the star-star relation \cite{Baxter:1986,Bazhanov:1992jqa} is precisely the
invariance of the 4d lens index
under Seiberg duality (see \cite{Spiridonov:2010em} for the case with $r=1$).
The star-star relation is stronger than the YBE and implies the latter.

We will identify the counterpart of YBE
to be a duality (called the Yang-Baxter duality)
between two 4d $\ensuremath{\mathcal{N}}=1$
quiver gauge theories, which can be obtained by composition/gauging of four
$\ensuremath{\mathcal{N}}=1$ Seiberg dualities.
Since this is a physical duality,
we expect that the discussion of integrability here
actually goes beyond the 4d lens indices discussed in this paper,
see comments in section \ref{sec.conclusion}.

The identification of integrable models and the 4d $\ensuremath{\mathcal{N}}=1$
gauge theories in the Gauge/YBE correspondence is summarized in Table \ref{tab.summary}.

\begin{table}[htbp]
\caption{Dictionary in the Gauge/YBE correspondence.}
\begin{center}
{\renewcommand\arraystretch{1.2}
\begin{tabular}{c|c}
integrable model & 4d $\ensuremath{\mathcal{N}}=1$ gauge theory \\
\hline
spin lattice & quiver diagram \\
rapidity line & global symmetry \\
rapidity (spectral) parameter & global symmetry fugacities \\
partition function & 4d lens index (on $S^1\times S^3/\ensuremath{\mathbb{Z}}_r$) \\
temperature-like parameters & fugacity $p, q$ \\
$U(1)$ spin variables & Wilson line along thermal $S^1$ \\
$\ensuremath{\mathbb{Z}}_r$ spin variables & Wilson line along the Hopf fiber \\
number of spin components & rank of a gauge group \\
self-interaction & $\ensuremath{\mathcal{N}}=1$ vector multiplet \\
nearest-neighbor interaction & $\ensuremath{\mathcal{N}}=1$ bifundamental chiral multiplet  \\
star-star relation & $\ensuremath{\mathcal{N}}=1$ Seiberg duality \\
R-matrix & theory $\mathcal{T}[R]$ (Figure \ref{fig.faceweight}) \\
composition of R-matrices & gauging \\
Yang-Baxter equation & Yang-Baxter duality \eqref{catYBE} \\
high-temperature expansion & dimensional reduction
\end{tabular}
}
\end{center}
\label{tab.summary}
\end{table}

\bigskip
The rest of this paper is organized as follows.
In section \ref{sec.model} we describe our statistical mechanical model,
whose gauge theory origin is explained in section \ref{sec.quiver}.
Further comments on the model are included in section
\ref{sec.comments},
and we conclude with summary and
comments on gauge theory implications in
section \ref{sec.conclusion}.
In Appendix \ref{sec.gamma}
we summarize the lens elliptic gamma function used in the paper.

We hope that the paper is readable both for
experts on supersymmetric gauge theories
and integrable models.
Experts on integrable models will
find the complete definition of the statistical
mechanical model in section \ref{sec.model},
without referring to the gauge theory explanation in
section \ref{sec.quiver}.
Experts on supersymmetric gauge theories
will benefit from section
\ref{sec.quiver}, which elucidates otherwise ad-hoc
definitions in section \ref{sec.model}.

\section{New Integrable Models}\label{sec.model}

In this section we spell out the definitions of our integrable lattice
spin models.
Our construction follows closely the $r=1$ case of
\cite{Bazhanov:2011mz} (see also \cite{Baxter:1997tn}).\footnote{However
some details are changed from \cite{Bazhanov:2011mz}
for the better match with supersymmetric gauge theories.
For example, the direction of
parallel rapidity lines (to be introduced momentarily)
are all the same in \cite{Bazhanov:2011mz}, whereas
they are alternating here. Our choice of rapidity parameters is
slightly more general than in \cite{Bazhanov:2011mz}.
}

\subsection{Definition of the Model}

\paragraph{Spin Lattice}

Let us consider a periodic spin lattice shown in Figure \ref{fig.lattice}.
We denote the set of vertices and edges of the lattice by $V$ and $E$.
Note that the edges in the lattice are oriented:
for an edge $e$ we denote its head (target) and tail (source) by $h(e), t(e)\in V$.
Note also that around each vertex the orientations of the arrows
alternate, and hence the number of incoming and outgoing arrows are the
same.
We color a vertex either black or white
depending on the orientations of the arrows around the vertex.
The lattice is then bipartite, {\it i.e.}\ an edge always connect
two vertices of different colors.\footnote{This bipartite graph is
{\it not} the bipartite graph discussed in
\cite{Yamazaki:2012cp,Terashima:2012cx};
it is the spin lattice (which will be identified with a
quiver diagram)
which is bipartite here.}
In this section we consider the square lattice,
see however comments on generalization towards the end of section \ref{sec.quiver}.

\begin{figure}[htbp]
\centering{\includegraphics[scale=0.3]{lattice.eps}}
\caption{The spin lattice of the statistical mechanical model.
The spins reside at the vertices of the lattice, and their
 nearest-neighbor interactions are represented by edges.
The interaction is chiral, which is represented by the orientation of
 the edges. The color of a vertex (black or white) is determined by the orientations
of arrows around it.}
\label{fig.lattice}
\end{figure}

\paragraph{Spin Variables}

We place spin variables
at each lattice site $v\in V$:
continuous $U(1)$ $N$-component spins
$z_{v,1}, \ldots, z_{v,N} \in e^{2\pi\sqrt{-1} \mathbb{R}} \, (|z_{v,i}|=1)$
 and discrete $N$-component spins $m_{v,1}, \ldots, m_{v,N} \in \ensuremath{\mathbb{Z}}_r$.
We moreover impose the condition that
\begin{equation}\begin{aligned}
\prod_i  z_{v,i}=1 \ , \qquad \sum_i m_{v,i}\equiv 0 \quad (\textrm{mod
} \, r) \ ,
\label{zprod}
\end{aligned}\end{equation}
and hence only $N-1$ of $z_{v,i}$ ($m_{v,i}$) are independent.
In the following we collectively denote the spins at a vertex $v$
by
$s_v=(z_v, m_v)$.

Note that we here consider classical (not quantum) statistical mechanical models,
and spin here just refers to dynamical variables of the theory over which we sum/integrate over.
On the boundary of the lattice we can impose periodic boundary conditions
or fixed boundary conditions. In the latter case the boundary spins will be non-dynamical
and the partition function will be a function of the values of the boundary spins.
For simplicity we in the following will mostly consider periodic
boundary conditions: the case with other boundary condition
is similar, but some part of the analysis requires careful analysis of boundary
effects.

\paragraph{Rapidity Parameters and Rapidity Lines}

For the purpose of discussing integrability let us introduce extra parameters to the
statistical model.
These parameters $R_e$ are associated for edges $e\in E$,
and we impose that condition that
\begin{equation}\begin{aligned}
\sum_{e: \textrm{ around a vertex } v}  (1-R_e)=2 \ , \quad
\sum_{e: \textrm{ around a face }}  R_e =2 \ .
\label{sumR}
\end{aligned}\end{equation}
It is these parameters $R_e$ which enter into the definition of the
partition function.

The conditions \eqref{sumR} do not have unique solutions.
However the ambiguity can be naturally
parametrized by the rapidity lines
\cite{Baxter:1978xr}\footnote{Rapidity
lines are called zig-zag paths in the context of
dimer models, which are the dual graphs of the spin lattice
discussed in this paper.}, which are drawn as
red line in Figure \ref{fig.latticewithlines}.
Rapidity lines contain the same data as the original
spin lattice --- we can recover the original lattice from the set of rapidity
lines, following the rules in Figure \ref{fig.latticerule}.

\begin{figure}[htbp]
\centering{\includegraphics[scale=0.3]{latticewithlines.eps}}
\caption{The rapidity lines for the spin lattice
of Figure \ref{fig.lattice}, shown in red arrows.
}
\label{fig.latticewithlines}
\end{figure}

\begin{figure}[htbp]
\centering{\includegraphics[scale=0.3]{latticerule.eps}}
\caption{We can reconstruct the spin lattice from the rapidity lines,
following the rule shown here.}
\label{fig.latticerule}
\end{figure}

Let us denote the set of rapidity lines by $I$,
and its element by $i\in I$. We associate a rapidity parameter
(spectral parameter)
$R_i$ to each rapidity line $i$.
Given a solution $R_e$ to \eqref{sumR}
we can modify $R_e$ to be\footnote{Since we have a
difference of two rapidities here
we have the so-called rapidity difference
property.}
\begin{equation}\begin{aligned}
R_e':=R_e+R_i-R_j \ ,
\label{difference}
\end{aligned}\end{equation}
when two rapidity lines $i, j$ pass through the edge $e$
as in Figure \ref{fig.rapidityrule}:
We can then easily verify that $R_e'$ also satisfy \eqref{sumR}.
Note that the overall shift of $R_i$
do not change the values of $R_e$, and
hence is irrelevant for the statistical mechanical model of this paper.

\begin{figure}[htbp]
\centering{\includegraphics[scale=0.3]{rapidityrule.eps}}
\caption{The rapidity parameter for an edge $e$ is given by the
 difference of $R_i$ and $R_j$ associated with two rapidity lines $i, j$.}
\label{fig.rapidityrule}
\end{figure}

\paragraph{Boltzmann Weight}

An edge $e\in E$, starting from a vertex $t(e)$ and ending on another $h(e)$,
represents the nearest-neighbor interaction
between the spins $s_{t(e)}$ and $s_{h(e)}$.
The corresponding Boltzmann weight is given by
\begin{equation}
\begin{split}
\ensuremath{\mathbb{W}}^e(R;s)&=\ensuremath{\mathbb{W}}^e_{R_e}(s_{t(e)}, s_{h(e)})\\
&=\mathcal{W}^e_0(s_{t(e)}, s_{h(e)})\, \prod_{1\le i,j\le N} \,
\Gamma_{r, \mmod{m_{t(e),i}-m_{h(e),j}}}\left(
(pq)^{\frac{R_e}{2}} \frac{z_{t(e),i}}{z_{h(e),j}}
;p,q\right) \ ,
\end{split}
\label{Wweight}
\end{equation}
where $\Gamma_{r,\mmod{m}}$ is the lens elliptic gamma function
and $\mmod{m}$ denotes $m$ modulo $r$ (see Appendix \ref{sec.gamma}).
The multiplicative factor $\mathcal{W}^e_0(s_{t(e)}, s_{h(e)})$ is given by
\begin{equation}
\begin{split}
 \mathcal{W}^e_0(s_{t(e)}, s_{h(e)})=&
 \prod_{1\le i,j \le N} \Bigg[ (p\, q)^{\frac{1}{4 r}\mmod{m_{t(e),i}-m_{h(e),j} } \mmod{-m_{t(e),i}+m_{h(e),j} }
(1-R_e) } \\
& \times
\, \left(\frac{p}{q} \right)^{-\frac{1}{12 r} \mmod{m_{t(e),i}-m_{h(e),j} } \mmod{-m_{t(e),i}+m_{h(e),j}}
(2 \mmod{m_{t(e),i}-m_{h(e),j} }-r)} \\
& \times  \left(\frac{z_{t(e),i}}{z_{h(e),j}} \right)^{-\frac{1}{2 r} \mmod{m_{t(e),i}-m_{h(e),j} } \mmod{-m_{t(e),i}+m_{h(e),j} }}
\Bigg]
\ .
\end{split}
\end{equation}
Note that this Boltzmann weight is chiral ({\it i.e.}\ not reflection
symmetric) $\ensuremath{\mathbb{W}}^e_{R_e}(s_{t(e)}, s_{h(e)})\ne \ensuremath{\mathbb{W}}^e_{R_e}(s_{h(e)},
s_{t(e)})$, except for the special case $N=2$.

We also include self-interactions among the
different components of the spin $s_v$ at the same vertex $v$.
The Boltzmann weight for the self-interaction at a vertex $v\in V$
is given by
\begin{equation}\begin{aligned}
\ensuremath{\mathbb{S}}^v(s)=\ensuremath{\mathbb{S}}^v(s_v)
=
\left( \prod_{a=0}^{r-1} \frac{1}{n_{v,a}!}
\right)
\mathcal{S}^v_0(s_v)
\left( \Gamma_{r, 0}(1;p,q) \right)^{-(N-1)}
\prod_{i \ne j }\, \Gamma_{r, \mmod{m_i-m_j}}\left(\frac{z_{v,i}}{z_{v,j}};p,q\right)^{-1}
\ ,
\label{Sweight}
\end{aligned}\end{equation}
where
we denoted the number of
$a \in \{0,1, \ldots, r-1\}$ in $\{m_{v,1}, \ldots, m_{v,N} \}\in
(\ensuremath{\mathbb{Z}}_r)^N$ to be $n_{v, a}$:
by definition we have $\sum_{a=0}^{k-1} n_{v,a}=N$.
The multiplicative factor $\mathcal{S}^v_0(s_v)$, which depends only on the
discrete spins $m_v$ out of $s_v$, is given by
\begin{equation}\begin{aligned}
\mathcal{S}^v_0(s_v)= (p \, q)^{-\frac{1}{4 r}\sum_{i\ne j} \mmod{m_{v,i}-m_{v,j} } \mmod{-m_{v,i}+m_{v,j} } } \ .
\end{aligned}\end{equation}
Note that the parameters $R_e$ do not appear in $\ensuremath{\mathbb{S}}^v$.

It is often possible to absorb this contribution into the
definition of $\ensuremath{\mathbb{W}}^e(s_{t(e)}, s_{h(e)})$, however
it will be more natural to keep this factor for the identification with gauge theories.

\paragraph{Partition Function}

The partition function is defined to be the statistical sum over all (continuous as well as
discrete) spin configurations $\{s_v \}_{v\in V}$:
\begin{equation}
\begin{split}
Z
&=\sum_{s_v\, (v\in V)} \,
\,
\left( \prod_{v\in V}  \ensuremath{\mathbb{S}}^v(s)\right)
\left( \prod_{e\in E} \ensuremath{\mathbb{W}}^e(R;s) \right) \\
&=\sum_{m_v\, (v\in V)} \int_{|z_{v,m}|=1}  \prod_{v\in V}
 \prod_{i=1}^{N-1}  \frac{d z_{v,i}}{2\pi \sqrt{-1} z_{v,i}}\
\left( \prod_{v\in V}  \ensuremath{\mathbb{S}}^v(z,m)\right)
\left( \prod_{e\in E} \ensuremath{\mathbb{W}}^e(R;z,m) \right) \ ,
\label{Ztotal}
\end{split}
\end{equation}
where in the integrand we write $z_{v,N}$ in terms of $z_{v,i}$ $(i=1, \ldots, N-1)$
in terms of \eqref{zprod}.

The total partition
function \eqref{Ztotal} depends on
two temperature-like parameters $p, q$ as well as rapidity parameters $R_i$.
All of these parameters will appear inside the Boltzmann weight of the
model. We can either regard $p, q$ as formal parameters, or impose $|p|, |q|<1$ for convergence of the Boltzmann
weights.

\paragraph{Reformulation as IRF Model}

The statistical mechanical model defined here
can be reformulated
either as a vertex model or
an IRF (interaction-round-a-face) model \cite{Bazhanov:2011mz}.
To avoid repetition we
in the following
use the language of IRF models.

To define an IRF model, we simply
integrate out spin variables
at white vertices,
while keeping those spins at black vertices.
The remaining lattice is then a larger lattice
colored green in
Figure \ref{fig.IRF},
and the Boltzmann weights are associated with the faces of the
new lattice.

Suppose that we have a face $F$,
with vertices $a, b, c, d$ as in Figure
\ref{fig.faceweight}. We denote the rapidity parameters of
four edges by $\alpha, \beta, \gamma, \delta$.
It follows from \eqref{sumR} that they satisfy the relation
$\alpha+\beta+\gamma+\delta=2$.
We take the Boltzmann weight for the face $F$
to be
\begin{equation}
\begin{split}
W(F)&=W\left(
\begin{array}{cc}
 d & c \\
 a & b \\
\end{array}
\right) \\
&=
\sqrt{\frac{\ensuremath{\mathbb{W}}_{2-\gamma-\delta}(d,c) \ensuremath{\mathbb{W}}_{2-\beta-\gamma}(b,c)}
{\ensuremath{\mathbb{W}}_{2-\alpha-\beta}(a,b)\ensuremath{\mathbb{W}}_{2-\alpha-\delta}(a,d)}}
\sqrt{\ensuremath{\mathbb{S}}^a \ensuremath{\mathbb{S}}^c}\, \sum_g
\ensuremath{\mathbb{S}}^g \, \ensuremath{\mathbb{W}}_{\alpha}(a,g) \ensuremath{\mathbb{W}}_{\beta}(g,b) \ensuremath{\mathbb{W}}_{\gamma}(c,g)
\ensuremath{\mathbb{W}}_{\delta}(g,d) \ ,
\end{split}
\label{Wdef}
\end{equation}
where we used the short-handed notation
$\ensuremath{\mathbb{W}}_{R_e}(t(e), h(e)):=\ensuremath{\mathbb{W}}^e_{R_e}(s_{t(e)}, s_{h(e)})$.

The partition function \eqref{Ztotal}
is can be re-written as
the partition function of this IRF model:
\begin{equation}\begin{aligned}
Z=\sum_{F\textrm{: face}} W(F) \ .
\label{ZF}
\end{aligned}\end{equation}
Note that we have also included
the Boltzmann weight for self-interactions $\ensuremath{\mathbb{S}}^v$
into the definition of the face weight \eqref{Wdef}.
In \eqref{Wdef} we have included a factor of
\begin{equation}\begin{aligned}
\sqrt{\frac{\ensuremath{\mathbb{W}}_{2-\gamma-\delta}(d,c) \ensuremath{\mathbb{W}}_{2-\beta-\gamma}(b,c)}
{\ensuremath{\mathbb{W}}_{2-\alpha-\beta}(a,b)\ensuremath{\mathbb{W}}_{2-\alpha-\delta}(a,d)}}
\label{extraW}
\end{aligned}\end{equation}
for later convenience.
This factor cancels out in the definition of the partition
function \eqref{ZF}, when we sum over faces of the IRF model.
Note that we have
$W_{R_e}(a,b)=W_{2-R_e}(b,a)^{-1}$
thanks to the relation \eqref{gammainv}.

\begin{figure}[htbp]
\centering{\includegraphics[scale=0.25]{IRF.eps}}
\caption{The reformulation as an IRF model. The spins at white
vertices are now integrated out, leaving only the black vertices
defined on the green lattice. The Boltzmann weight $W$ is
then associated with the faces of the new lattice.
Alternative possibility is to exchange the role of white and black
 vertices,
leading to another IRF model defined on the dotted lattice.
The two IRF models nevertheless have the same partition function.
}
\label{fig.IRF}
\end{figure}

\begin{figure}[htbp]
\centering{\includegraphics[scale=0.3]{faceweight.eps}}
\caption{The face weight for the IRF model is determined
from the Boltzmann weights for four edges, as well
as those of four vertices, of the original spin lattice.
We have denoted the four black vertices by $a, b, c, d$,
and the weights for the edges connected by them by
$\alpha, \beta, \gamma, \delta$. The first of the constraint equations \eqref{sumR}
in this case means $\alpha+\beta+\gamma+\delta=2$.
}
\label{fig.faceweight}
\end{figure}

\subsection{Integrability}\label{sec.lattice}

We now claim that the model discussed above
is integrable.
More precisely the model satisfies the
star-star relation and consequently the YBE.
This is the main result of the paper.

\paragraph{YBE}

The integrability is defined by the famous Yang-Baxter equation.
For the IRF model, the YBE is written as \cite{Baxter:1982zz}
\begin{equation}
\begin{split}
\sum_g \matW{\alpha}{\beta}{\eta}{\zeta}{a}{b}{g}{f} \matW{\zeta}{\eta}{\delta}{\epsilon}{f}{g}{d}{e}
\matW{\eta}{\beta}{\gamma}{\delta}{g}{b}{c}{d}
=\sum_g \matW{\zeta}{\alpha}{\eta}{\epsilon}{f}{a}{g}{e}
 \matW{\alpha}{\beta}{\gamma}{\eta}{a}{b}{c}{g}
\matW{\eta}{\gamma}{\delta}{\epsilon}{g}{c}{d}{e} \
 ,
\end{split}
\end{equation}
which can be graphically represented as in
Figure \ref{fig.YBE}(a).
The YBE implies the commutativity of the row-to-row transfer
matrices, leading to an infinite number of conserved charges \cite{Baxter:1982zz}.

The YBE can also be written in terms of the
original Boltzmann weights $\ensuremath{\mathbb{S}}^v, \ensuremath{\mathbb{W}}^e$
(recall the relation \eqref{Wdef}),
which gives the relation in Figure \ref{fig.YBE}(b).
Most of the contributions from \eqref{extraW}
cancel out, however some of them remain and appear in
Figure \ref{fig.YBE}(b) as extra arrows
connecting two black vertices.

\begin{figure}[htbp]
\centering{\includegraphics[scale=0.25]{YBE.eps}}
\caption{YBE for the IRF model (a) and for the rapidity lines (b).}
\label{fig.YBE}
\end{figure}

\paragraph{Star-Star Relation}

Instead of solving the YBE directly,
we solve the the star-star (also called star-to-reversed-star) relation
\cite{Baxter:1986,Bazhanov:1992jqa}:\footnote{
In the definition of the IRF model,
we can choose to integrate out black
vertices instead of white vertices.
This leads to a different IRF model with a face weight
\begin{equation}\begin{aligned}
\overline{W}(F)=\overline{W}
\left(
\begin{array}{cc}
 d & c \\
 a & b \\
\end{array}
\right)
=
\sqrt{\frac{\ensuremath{\mathbb{W}}_{2-\alpha-\beta}(a,b)\ensuremath{\mathbb{W}}_{2-\alpha-\delta}(a,d)}{\ensuremath{\mathbb{W}}_{2-\gamma-\delta}(d,c) \ensuremath{\mathbb{W}}_{2-\beta-\gamma}(b,c)}
}
\sqrt{\ensuremath{\mathbb{S}}^a \ensuremath{\mathbb{S}}^c}\, \sum_g
\ensuremath{\mathbb{S}}^g \,
\ensuremath{\mathbb{W}}_{\gamma}(g,a) \ensuremath{\mathbb{W}}_{\delta}(b,g) \ensuremath{\mathbb{W}}_{\alpha}(g,c)
\ensuremath{\mathbb{W}}_{\beta}(d,g) \ .
\end{aligned}\end{equation}
The partition function is again given by \eqref{ZF},
where this time the sum is over the faces of another lattice with
white vertices only.
The star-star relation \eqref{starstar}
ensures that the two Boltzmann weights
$W$ and $\overline{W}$
are the same.
}
\begin{equation}
\begin{split}
&\ensuremath{\mathbb{W}}_{2-\gamma-\delta}(d,c) \ensuremath{\mathbb{W}}_{2-\beta-\gamma}(b,c)
\sum_g \ensuremath{\mathbb{S}}^g\, \ensuremath{\mathbb{W}}_{\alpha}(a,g) \ensuremath{\mathbb{W}}_{\beta}(g,b)
\ensuremath{\mathbb{W}}_{\gamma}(c,g) \ensuremath{\mathbb{W}}_{\delta}(g,d) \\
&\quad= \ensuremath{\mathbb{W}}_{2-\alpha-\beta}(a,b) \ensuremath{\mathbb{W}}_{2-\alpha-\delta}(a,d)
\sum_g \ensuremath{\mathbb{S}}^g\, \ensuremath{\mathbb{W}}_{\gamma}(g,a) \ensuremath{\mathbb{W}}_{\delta}(b,g)
\ensuremath{\mathbb{W}}_{\alpha}(g,c) \ensuremath{\mathbb{W}}_{\beta}(d,g) \ .
\end{split}
\label{starstar}
\end{equation}
The graphical representations of this relation is given in Figure \ref{fig.starstar}.\footnote{This
was called the double Yang-Baxter move in \cite{Yamazaki:2012cp};
it was called the ``double'' since the Yang-Baxter equation is
often associated with a crossing of three lines, and the star-star relation involves
two such moves and four lines.
Note however the YBE of our statistical mechanical model can be represented
as re-shuffling to three sets of parallel lines (Figure \ref{fig.YBE}),
totally six lines.
In integrable model language our statistical mechanical model satisfies the star-star relation,
but not the star-triangle relation.
}
We can directly verify that the YBE follows from the star-star
relation, see Figure \ref{fig.YBEfromStar}.

\begin{figure}[htbp]
\centering{\includegraphics[scale=0.3]{starstar.eps}}
\caption{The graphical representation of the star-star relation.}
\label{fig.starstar}
\end{figure}

\begin{figure}[htbp]
\centering{\includegraphics[scale=0.25]{YBEfromStar.eps}}
\caption{YBE (Figure \ref{fig.YBE}) follows from a sequence of star-star
 relations (Figure \ref{fig.starstar}) \cite{Baxter:1997tn}:
in each step the star-star relation is taken around the blue vertex.}
\label{fig.YBEfromStar}
\end{figure}

It would be desirable to mathematically prove the star-star relation.
For the case with $r=1$ and general $N$, the star-star relation follows from the results of
\cite{RainsTransf}, as first pointed out by \cite{Dolan:2008qi}
(see also \cite{Romelsberger:2007ec} for earlier work).\footnote{For the $r=1, N=2$
case (where the quiver gauge theory is non-chiral),
the star-star relation follows from the
simpler relation, the star-triangle relation,
which can be proven \cite{Spiridonov:2010em} from the
elliptic beta integral of \cite{SpiridonovBeta}.}
In this paper the star-star relation will be derived from the
considerations in
supersymmetric gauge theory.
This is be the topic of the next section.

\section{Gauge Theory Origin}\label{sec.quiver}

Let us now comment
on the gauge-theoretic origin of our
statistical mechanical model.
Interestingly, each of the ingredients introduced in the previous section
has a precise counterpart in
4d $\ensuremath{\mathcal{N}}=1$ supersymmetric gauge theories
(see Table \ref{tab.summary} in section \ref{sec.intro}).
The gauge theory moreover gives generalization of the
statistical mechanical model to spin lattices
more general than the square lattices of the previous section.
The discussion here will be brief since
most of the ingredients here have been
discussed in detail in
\cite{Yamazaki:2012cp,Terashima:2012cx,Xie:2012mr}.

\subsection{Star-Star Relation as Seiberg Duality}

The first step is to define a 4d $\ensuremath{\mathcal{N}}=1$ theory associated with the
spin lattice.

The idea is to regard the spin lattice as
the quiver diagram.  Namely we associate an
$SU(N)$ gauge group to a vertex $v\in V$,
a bifundamental field to an edge $e\in E$,
and a superpotential term (a trace of the product of bifundamentals)
to a face of the spin lattice.
Note that the orientation of the graph corresponds to the chirality of bifundamental
multiplets.

It is believed in the literature that the resulting model flows in the IR
to a non-trivial strongly-coupled fixed point. We assume in this paper that this is the
case.\footnote{For periodic boundary
conditions there is a strong support for this
from the AdS/CFT correspondence (see \cite{Kennaway:2007tq,Yamazaki:2008bt}
and references therein).
See also \cite{Heckman:2012jh}
for the evidence
of non-trivial IR fixed points for the theories of \cite{Xie:2012mr}.
}
Since the theory has many $U(1)$ global symmetries
the $U(1)$ R-symmetry in the IR can be a mixture of UV
$U(1)$ R-symmetry with other global $U(1)$ symmetries,
leading to an ambiguity of
the choice of the R-charge $R_e$
for a bifundamental at an edge $e\in E$.
The ambiguities of the IR R-charges coming from the mixing with
global symmetries
can be parametrized as in \eqref{difference} \cite{Hanany:2005ss}\footnote{In
\cite{Hanany:2005ss} the rapidity parameter was identified with the
values of the IR R-charge in the superconformal algebra
determined from $a$-maximization.
However for our purposes we do not adopt this interpretation,
and rather it is crucial to use allow
the rapidity values to be away from their IR values.
}, where $R_i$ represents the global $U(1)$ symmetry
associated with a rapidity line.
The conditions
\eqref{sumR} on the R-charges
represent the vanishing of the beta-functions
for the gauge couplings and superpotential couplings, respectively.

Now
the crucial observation (which goes back to \cite{Spiridonov:2010em}) is that the star-star relation \eqref{starstar}, when translated into the language of
quiver gauge theories, is precisely the 4d $\ensuremath{\mathcal{N}}=1$ Seiberg duality \cite{Seiberg:1994pq}--- in Figure \ref{fig.starstar}
we have Seiberg duality for $N_f=2 N$, where the global symmetry
$SU(2N)$ is broken to $SU(N)\times SU(N)$ by a
superpotential.\footnote{Mathematically the change of
quiver in Figure \ref{fig.starstar}
is an example of a quiver mutation.
}
The invariance of the partition function under star-star
relation is then translated into the invariance of
some quantity under 4d Seiberg duality.

There is one nice quantity which is indeed invariant under the
Seiberg duality: the 4d lens index, the twisted partition function on
$S^1\times S^3/\ensuremath{\mathbb{Z}}_n$. Here $S^3/\ensuremath{\mathbb{Z}}_n$ is the lens space,
where $\ensuremath{\mathbb{Z}}_n$ acts on the Hopf fiber of $S^3$.
Let us now discuss this quantity.

\subsection{Partition Function from 4d Lens Index}

The 4d lens index \cite{Benini:2011nc} is defined by
\begin{equation}\begin{aligned}
\ensuremath{\mathcal{I}}(p,q,u_i)=\textrm{Tr}_{\, \ensuremath{\mathcal{H}}_{S^3/\ensuremath{\mathbb{Z}}_r}} \left[
(-1)^F p^{\frac{E+j_2}{3}+j_1} q^{\frac{E+j_2}{3}-j_1}
\right] \ ,
\label{Idef}
\end{aligned}\end{equation}
Here $\ensuremath{\mathcal{H}}_{S^3/\ensuremath{\mathbb{Z}}_r}$ is the Hilbert space on
$\ensuremath{\mathcal{H}}_{S^3/\ensuremath{\mathbb{Z}}_r}$,
$j_1, j_2$ are the generators of $U(1)\times U(1)\in SU(2)\times
SU(2)\sim SO(4)$ and $E$ the generator of translation along $S^1$.
The fugacities for the global symmetries, which can be included
in the definition of \eqref{Idef}, is traded for
the choice of R-charges.
In other words the fugacity for the $i$-th global symmetry is
given by $(pq)^{R_i}$.

The lens index is one of the most powerful quantities
known in the literature for
quantitative analysis of
4d $\ensuremath{\mathcal{N}}\ge 1$ supersymmetric gauge theories.
It is more powerful than
their $r=1$ counterparts (see \cite{Razamat:2013jxa,RazamatWillett}
for examples), and
also reproduce many partition functions
in lower dimensions, such as the
3d lens space partition function, the
3d superconformal index and the 2d sphere partition function
(\cite{Yamazaki:2013fva} and references therein).

Since the 4d lens index is an index, it is invariant under continuous deformations and
is invariant under the RG flow, as long as we take into account the change of the R-charge
under the RG flow ({\it cf.}\ \cite{Romelsberger:2007ec}).\footnote{In principle the same logic
applies to any 4d index invariant under Seiberg duality,
however at present the 4d lens index
is most general index known in the literature.}
This means that
the index can be evaluated in the free field limit.
The resulting expression coincides with
expression found in section \ref{sec.model}.
For example,
the spin variables $s_v=(z_v, m_v)$ are the Wilson lines
for the $SU(N)$ gauge group at vertex $v\in V$.
Since $\pi_1(S^1\times S^3/\ensuremath{\mathbb{Z}}_r)=\ensuremath{\mathbb{Z}}\times \ensuremath{\mathbb{Z}}_r$,
we have continuous Wilson lines
$z_v=(z_{v,1}, \ldots, z_{v,N})\in U(1)^N$,
(with $\prod_i z_{v,i}=1$)
along the $S^1$,
or discrete Wilson lines
$m_v=(m_{v,1}, \ldots, m_{v, N})\in (\ensuremath{\mathbb{Z}}_r)^{N}$
(with $\sum_i m_{v,i}=0$ modulo $r$)
in the Hopf fiber direction of the manifold.
This explains our choice of spin variables.

The Boltzmann weights $\ensuremath{\mathbb{S}}^v, \ensuremath{\mathbb{W}}^e$
are the 1-loop determinants for
vector and chiral multiplets
(including the zero-point contributions
$\mathcal{S}^v_0, \mathcal{W}^e_0$),
whose explicit expression can
be found in \cite{Benini:2011nc,Yamazaki:2013fva}.\footnote{For 4d $\ensuremath{\mathcal{N}}=1$ theories
there was a sign error in version 2 of \cite{Benini:2011nc}, which is corrected in the current paper.
The author would like to thank S.S. Razamat and B. Willet for pointing this out.
He would also like to thank S.S. Razamat for discussion on zero-point contributions
of the 4d lens index.
}
The chiral multiplet 1-loop determinant $\ensuremath{\mathbb{W}}^e$
corresponds to a nearest-neighbor interaction on the spin lattice
since a bifundamental is charged only with respect to
two $SU(N)$ gauge groups.

The star-star relation represents the Seiberg duality.\footnote{See
\cite{RazamatWillett} for checks of Seiberg duality for 4d lens indices.
}
The factors
\begin{equation}\begin{aligned}
\ensuremath{\mathbb{W}}_{2-\gamma-\delta}(d,c), \quad \ensuremath{\mathbb{W}}_{2-\beta-\gamma}(b,c),\quad
\ensuremath{\mathbb{W}}_{2-\alpha-\beta}(a,b), \quad \ensuremath{\mathbb{W}}_{2-\alpha-\delta}(a,d)
\label{meson}
\end{aligned}\end{equation}
in \eqref{starstar}, which can be traced back to \eqref{extraW},
represents the contributions from the mesons in the Seiberg duality.
Since a meson $M$ couples to a quark $q$ and anti-quark $\tilde{q}$
by a superpotential term $W=\tilde{q}Mq$,
and hence their R-charges should sum up to $2$.
This is represented by the parametrizations of the rapidity parameters
in \eqref{meson} (in particular the numbers ``2'' there).

The statistical model of section \ref{sec.lattice}
uses the square lattice (with periodic boundary conditions).
However the gauge theory considerations apply to more general lattices
determined from the rapidity lines,
as long as the quiver defined from them is connected
(more precisely this is the
the admissibility conditions of \cite{Ueda:2006jn,Yamazaki:2012cp}).
The vertex of the lattice can be
of $2n$-valent ($n\ge 2$).
The discussion of the star-star relation or equivalently the Seiberg
duality
works in exactly the same manner.

\section{Further Comments}\label{sec.comments}

Several comments are now in order.

\paragraph{Specialization and Limits}

Since the star-star relation is an equality,
the limits and specializations of the model is still guaranteed to be
integrable.
We can for example take the fugacities $p,q$
to be at root of unity, or regard them as temperatures
and take the high-temperature limit.

We can for example take $r=1$.
The discrete spins $m_v$
are then frozen to be trivial values $m_v=0$, and
the only non-trivial
degrees of freedom are continuous $U(1)$
spins $z_{v,1}, \ldots, z_{v,N}$.
As discussed already in introduction, this
coincides with the Bazhanov-Sergeev model \cite{Bazhanov:2011mz}
(up to the normalization constant of the partition function
which is not essential for the discussion of integrability here).
This model is known to reproduce a number of
known integrable models \cite{Bazhanov:2007mh,Bazhanov:2010kz,Bazhanov:2011mz}, such as
the Faddeev-Volkov model \cite{Volkov:1992uv,FaddeevVolkovAbelian}\footnote{The reduction
from the Bazhanov-Sergeev model to
Faddeev-Volkov model
corresponds to
the dimensional reduction and Higgsing of
the 4d $\ensuremath{\mathcal{N}}=1$ theory down to 3d $\ensuremath{\mathcal{N}}=2$ theory.
In this process 2d spin system lifts to the geometry of
a hyperbolic 3-manifold \cite{Yamazaki:2012cp,Terashima:2012cx}.
}, the chiral Potts model \cite{AuYang:1987zc,Baxter:1987eq} and the
 Kashiwara-Miwa model \cite{Kashiwara:1986tu}.
The last two models are different generalizations of the
Fateev-Zamolodchikov model \cite{Fateev:1982wi}.

Interestingly,
many of these models have discontinuous spins:
the values of the continuous spins in the Bazhanov-Sergeev model
are frozen to discrete values--- this happens, for example, when
we choose the temperature-like parameters $p,q$
to be at root of unity.
It is interesting to ask whether the inverse procedure is possible,
namely whether we can construct a model with all spins continuous
such that our model is reproduced as a high-temperature limit.
Such a model, if it exists, should be more powerful than the model discussed in
this paper and could have some gauge theory implications, for example
as an index of some 5d gauge theory.

\paragraph{Further Generalizations}

There are several natural generalizations of our model.
For example, we can take the
rank of gauge groups
(and hence the number of components of spins)
to be $v$-dependent: the gauge group at $v$ is $SU(N_v)$.
It is straightforward to write down the corresponding statistical model.
However, not all the possible values of $N_v$ are allowed
since we have the ensure the absence of anomalies
and the existence of the non-trivial IR fixed point.
Moreover the ranks of the gauge groups will change under the Seiberg
duality, and for this reason the corresponding YBE is more general than the one
typically considered in integrable models.
The same subtlety happens when we consider
$S^1\times S^2$ index of 3d $\ensuremath{\mathcal{N}}=2$ quiver Chern-Simons-matter theories
and the Giveon-Kutasov dualities \cite{Giveon:2008zn} among them.

\paragraph{Comparison with Gauge/Bethe Correspondence}

The Gauge/YBE correspondence \eqref{GaugeYBE} discussed in this paper
is reminiscent of a
similar correspondence in the literature,
the Gauge/Bethe correspondence \cite{Nekrasov:2009uh,Nekrasov:2009ui}:
\begin{equation}\begin{aligned}
\ensuremath{\mathcal{W}}_{\rm eff}(\sigma)=Y(\sigma) \ .
\label{GaugeBethe}
\end{aligned}\end{equation}
On the left hand side of \eqref{GaugeBethe}
we have the effective twisted
super potential $\ensuremath{\mathcal{W}}_{\rm eff}(\sigma)$
of a 2d $\ensuremath{\mathcal{N}}=(2,2)$ theory as a function of
the scalar in the twisted multiplet $\Sigma=\overline{D}_+ D_- V$.
On the right hand side we have the Yang-Yang function \cite{Yang:1968rm}
$Y(\sigma)$. The derivative of this function gives the Bethe Ansatz equation,
which coincides with the vacuum equation of the
2d $\ensuremath{\mathcal{N}}=(2,2)$ theory:
\begin{equation}\begin{aligned}
\exp\left(\frac{\partial \ensuremath{\mathcal{W}}_{\rm eff}}{\partial \sigma} \right)=
\exp\left(\frac{\partial Y}{\partial \sigma}\right)=1\ .
\end{aligned}\end{equation}

While the two relations \eqref{GaugeYBE}, \eqref{GaugeBethe}
are similar, there are crucial differences between the
two.\footnote{Despite
the differences it would be interesting to ask if the two relations
are related more directly, at least for a special class of theories.
In fact, if we choose $r=1$
the dimensional reduction of the 4d index
in the Gauge/YBE correspondence
gives the 3d ellipsoid partition function $S^3_b$,
and then the 2d effective twisted superpotential in the
$b\to 0$ limit. This reduction preserves integrability
of the associated statistical mechanical model.
Since we have 2d $\ensuremath{\mathcal{N}}=(2,2)$ theory,
we could then consider the
corresponding integrable model of the 2d theory through the Gauge/Bethe correspondence, if it exists at all.}
One crucial difference is that
in Gauge/YBE we have the partition function
of an integrable model,
whereas in Gauge/Bethe we have
the Yang-Yang function.
This is a huge difference when we want to
find {\it new} integrable models from the correspondence.

In Gauge/Bethe correspondence
it is in general
rather hard to identify the corresponding integrable model for
a given 2d $(2,2)$ theory---there is no general algorithm
to recover the R-matrices from
the Bethe Ansatz equation.
In fact, it is not known in general whether such an integrable model
really exists for a given 2d $\ensuremath{\mathcal{N}}=(2,2)$ theory.\footnote{The
exception
is the case where the 2d $(2,2)$ theory comes from 4d $\ensuremath{\mathcal{N}}=2$
class $\mathcal{S}$ theories
on the Omega-background \cite{Nekrasov:2002qd} with
equivariant parameters $\epsilon_1=0, \epsilon_2\ne 0$ \cite{Nekrasov:2009rc}:
the corresponding integrable model in this case
is Hitchin integrable model.
However, this is not really a new integrable model
found from the correspondence,
and moreover the story does not generalize to more general 2d $\ensuremath{\mathcal{N}}=(2,2)$ theories.}

This should be contrasted with the Gauge/YBE correspondence
discussed in this paper.
There the gauge theory answer directly gives
precisely the R-matrix and the partition function of
the integrable model, and the
the reason {\it why} the model is integrable
has a direct explanation from the 4d Seiberg duality.
This direct relation was the reason that in this paper
we could construct new integrable models
from the correspondence.

\paragraph{Underlying Algebra}

One of the systematic methods to construct solutions to YBE is to
define the R-matrix as an intertwiner of the tensor product of representations.
The relevant algebra for the $r=1$ case is the Sklyanin algebra \cite{Sklyanin:1982tf}
and their generalizations \cite{CherednikGeneralized}
(see e.g.\ \cite{Zabrodin:2010qm}), which are 2-parameter
deformations of the universal enveloping algebra $U(\mathfrak{gl}_N)$.
Such an algebra acts on the operators of {\it a class of} supersymmetric gauge theories,
and plays the role similar to that of the Yangian in the Gauge/Bethe correspondence \cite{Nekrasov:2009uh,Nekrasov:2009ui}.

\paragraph{Knot Invariants}

It would be interesting to study the knots invariants
associated with our solutions to the YBE ({\it cf.}\ \cite{Wadati:1989ud}),
see \cite{Wu2} for related discussion for the chiral Potts model.

\section{Concluding Remarks}\label{sec.conclusion}

We proposed a 2d statistical
spin lattice model with
nearest-neighbor and self interactions.
At each vertex we have a hybrid of continuous spins and discrete spins,
each of which has $N$ components.
The Boltzmann weights are given in \eqref{Wweight}, \eqref{Sweight}
and are written in terms of lens elliptic gamma functions \eqref{egamma1}.
The model can be reformulated either as a vertex model
or an IRF model.

The model is integrable (satisfy the star-star relation and hence the YBE)
if we assume
the invariance of the 4d lens index under the Seiberg duality.
Alternatively, the mathematical proof of integrability
can be thought of as another non-trivial test of
the existence of non-trivial IR fixed points
and the 4d Seiberg duality.

It is important to mathematically prove the
integrability of the model. We can also study
various properties of the model, such as free energies and correlation
functions, and their gauge theory
interpretations.\footnote{Thermodynamics limit of the model
could have direct geometric
interpretation ({\it cf.}\ \cite[section 5]{Yamazaki:2012cp}),
perhaps along the lines of \cite{KenyonOkounkovSheffield,Ooguri:2009ri}.
}

\bigskip

In this paper
we have concentrated on the application of supersymmetric
gauge theories to integrable models.
We hope that our results are of interest to
experts on integrable models.
While this is interesting in its own right,
experts on supersymmetric gauge theorists might
be interested in a different question,
namely if our results sheds any new light on
our understanding of
4d $\ensuremath{\mathcal{N}}=1$ supersymmetric gauge theories.

In this respect, one fundamental question is whether
the integrable structure discussed in this paper goes
beyond the 4d lens index, and extends to the
full physics of the IR fixed points.

We do not have the complete answer yet, however
let us here point out that
we have already identified dualities of quiver gauge theories
corresponding to YBE: it is the sequence of four Seiberg dualities,
and the quiver diagrams for the two theories
are shown in Figure \ref{fig.YBE}.
This duality among supersymmetric gauge theories
requires several $SU(N)$
gauge groups, which is probably the reason why
it has not been paid much attention so far.

Recall that the R-matrix in integrable models
is the linear map $R: V\otimes V \to V\otimes V$.
This corresponds to a 4d $\ensuremath{\mathcal{N}}=1$ theory
$\mathcal{T}[R]$, whose quiver diagram is a simply a ``star''
of Figure \ref{fig.faceweight}---
if we un-gauge the gauge groups at the four black vertices,
the theory is has global symmetry $SU(N)^4$,
and the $SU(N)$ global symmetry plays the role
similar to that of $V$.
The counterpart of taking a product of R-matrices
is to concatenate two quiver diagrams,
by re-gauging the diagonal part of global $SU(N)\times SU(N)$
symmetries of the theories associated with the quivers.
Finally, the counterpart of YBE is
the following duality between two 4d $\ensuremath{\mathcal{N}}=1$
theories:
\begin{equation}\begin{aligned}
\mathcal{T}[R_{12}] \circ \mathcal{T}[R_{13}] \circ \mathcal{T}[R_{23}]
\quad
\textrm{is dual to}
\quad \mathcal{T}[R_{23}] \circ \mathcal{T}[R_{13}] \circ
\mathcal{T}[R_{12}] \ ,
\label{catYBE}
\end{aligned}\end{equation}
where the symbol $\circ$
represents the gauging the diagonal subgroup of
product global symmetries of two theories.\footnote{More precisely
we need to include two extra bifundamentals chiral
multiplets (mesons)
to both sides
of \eqref{catYBE} (Figure \ref{fig.YBE}).
}
The duality \eqref{catYBE}, which we call the {\it Yang-Baxter duality},
is more powerful than the equivalence of 4d lens indices.
For example, the physical Hilbert spaces of the two theories should
coincide in the IR, which can be thought of
as a version of categorification of the Gauge/YBE correspondence.

The identification of the YBE with the
duality \eqref{catYBE}
makes it clear that for the discussion of
integrability we need to consider
not a single theory, but a class of them related by
gauging/un-gauging procedures:
in the spirit of \cite{Yamazaki:2013xva},
what is relevant here is ``integrability in theory space''.

It would be interesting to analyze the duality \eqref{catYBE}
further. The fact that the YBE is more general than the
star-star relation in integrable model suggests that
there should be generalizations of theory $\mathcal{T}[R]$
and dualities among them,
which cannot be decomposed into Seiberg dualities.
Such dualities will give an ultimate reason
{\it why} gauge-theoretic quantities have integrable structures.

\section*{Acknowledgments}

The author would like to take this opportunity to
thank M. Jimbo for his excellent textbook ``quantum groups and Yang-Baxter equations'' (Springer, 1990, in Japanese).
The author has been mesmerized by the beauty of integrable models
ever since he read the textbook a decade ago.
The author would also like thank S. S. Razamat and B. Willett
for stimulating discussion on the 4d lens index and for informing the author of
their forthcoming paper \cite{RazamatWillett}.
He would like to thank Yukawa Institute for Theoretical Physics, Kyoto University
(YKIS 2012) for hospitality where part of this work was performed.
The results of this paper was presented at the pre-strings meeting
``Exact Results in String/M-theory'' (KIAS, June 2013) and
the author would like to thank the audience for feedback.
This research was supported by World Premier
International Research Center Initiative (WPI Initiative), MEXT, Japan.