\section{Introduction}
\label{sec:1}
Ideally, a non-perturbative formulation of a field theory with some symmetries
should provide, not only the definition of correlation functions, but also the
definition of renormalized Noether currents that generate correctly-normalized
symmetry transformations on renormalized fields. This is expressed by
renormalized Ward--Takahashi (WT) relations and, when these relations hold, one
may say that the symmetries are \emph{really\/} realized in quantum field
theory.

Quite often, however, the regularization procedure breaks the preferred
symmetries and, for this reason, it is generally very difficult to conclude
even the existence of such renormalized Noether currents (especially when the
symmetry transformations are non-linear). Even if such Noether currents are
assumed to exist, the explicit construction can be very cumbersome. In the
lattice formulation of supersymmetric theories, one encounters such a
situation, because the infinitesimal translation (that is a part of the SUSY
algebra) is broken by the spacetime lattice, implying also the breaking of
SUSY. Here, almost all fundamental symmetries that define the theory are
broken by the regularization.

Having the above general remark in mind, in the present paper, we study the
construction of Noether currents in the lattice formulation of the
four-dimensional $\mathcal{N}=1$ supersymmetric Yang--Mills theory (4D
$\mathcal{N}=1$ SYM)~\cite{Curci:1986sm} (see also Ref.~\cite{Kaplan:1983sk}
for an earlier consideration and Ref.~\cite{Montvay:2001aj} for a very readable
review). The classical Euclidean action of 4D $\mathcal{N}=1$ SYM is given
by\footnote{We basically follow the notation of~Ref.~\cite{Suzuki:2012pc}:
The sum over repeated indices is understood. Vector indices $\mu$, $\nu$,
\dots, run over $0$, $1$, $2$, $3$. $\epsilon_{\mu\nu\rho\sigma}$~denotes the
totally anti-symmetric tensor and~$\epsilon_{0123}=-1$. All gamma matrices are
hermitian and obey $\{\gamma_\mu,\gamma_\nu\}=2\delta_{\mu\nu}$. We define
$\gamma_5\equiv-\gamma_0\gamma_1\gamma_2\gamma_3$
and~$\sigma_{\mu\nu}\equiv[\gamma_\mu,\gamma_\nu]/2$. The charge conjugation
matrix~$C$ satisfies, $C^{-1}\gamma_\mu C=-\gamma_\mu^T$,
$C^{-1}\sigma_{\mu\nu}C=-\sigma_{\mu\nu}^T$, $C^{-1}\gamma_5C=\gamma_5^T$
and~$C^T=-C$. The generator of the gauge group $SU(N_c)$, $T^a$, is normalized
as $\tr(T^aT^b)=(1/2)\delta^{ab}$. $F_{\mu\nu}(x)$ is the field strength
$F_{\mu\nu}(x)\equiv\partial_\mu A_\nu(x)-\partial_\nu A_\mu(x)%
+ig[A_\mu(x),A_\nu(x)]$, where $g$~is the bare gauge coupling constant, and
$\Slash{D}\equiv\gamma_\mu D_\mu$, where
$D_\mu\equiv\partial_\mu+ig[A_\mu(x),]$ is the covariant derivative in the
adjoint representation. For the lattice theory, $x$, $y$, $z$, \dots\ denote
lattice points and $a$~is the lattice spacing; $\Hat{\mu}$ is the unit vector
in the $\mu$-direction. The link variable~$U_\mu(x)$ and the gauge
potential~$A_\mu(x)$ are related as
\begin{equation}
   U_\mu(x)=e^{iagA_\mu(x)}.
\label{eq:(1.1)}
\end{equation}
The forward and backward difference operators respectively are defined by
\begin{equation}
   \partial_\mu f(x)
   \equiv\frac{1}{a}\left[f(x+a\Hat{\mu})-f(x)\right],\qquad
   \partial_\mu^*f(x)
   \equiv\frac{1}{a}\left[f(x)-f(x-a\Hat{\mu})\right],
\label{eq:(1.2)}
\end{equation}
and the symmetric difference operator $\partial_\mu^S$ is defined by
\begin{equation}
   \partial_\mu^S\equiv\frac{1}{2}(\partial_\mu+\partial_\mu^*).
\label{eq:(1.3)}
\end{equation}}
\begin{equation}
   S_{\text{classical}}
   =\int d^4x\,\left\{
   \frac{1}{2}\tr\left[F_{\mu\nu}(x)F_{\mu\nu}(x)\right]
   +\tr\left[\Bar{\psi}(x)\Slash{D}\psi(x)\right]\right\},
\label{eq:(1.4)}
\end{equation}
where the adjoint fermion (gluino)~$\psi(x)$ on the Euclidean space is subject
of the constraint
\begin{equation}
   \Bar{\psi}(x)=\psi^T(x)(-C^{-1}),
\label{eq:(1.5)}
\end{equation}
to express the degrees of freedom of a Majorana fermion in the 4D Minkowski
space. The global SUSY~$\Bar{\delta}_\xi$ in the classical 4D $\mathcal{N}=1$
SYM is
\begin{align}
   \Bar{\delta}_\xi A_\mu(x)&=\Bar{\xi}\gamma_\mu\psi(x),&&
\notag\\
   \Bar{\delta}_\xi\psi(x)&=-\frac{1}{2}\sigma_{\mu\nu}\xi F_{\mu\nu}(x),&
   \Bar{\delta}_\xi\Bar{\psi}(x)
   &=\frac{1}{2}\Bar{\xi}\sigma_{\mu\nu}F_{\mu\nu}(x),
\label{eq:(1.6)}
\end{align}
where the Grassmann-odd constant spinor~$\xi$ obeys $\Bar{\xi}=\xi^T(-C^{-1})$.
Without containing the scalar field, 4D $\mathcal{N}=1$ SYM is the simplest
supersymmetric theory in four dimensions and we should understand the issue
raised at the beginning of this paper first in this example, before tackling
more complicated supersymmetric theories.

The most interesting result we will obtain below is a very explicit form of an
energy-momentum tensor on the lattice $\mathcal{T}_{\mu\nu}(x)$, given
by~Eqs.~\eqref{eq:(3.3)} and~\eqref{eq:(3.4)}, which is \emph{conserved\/} in
the quantum continuum limit. Our way of construction of this energy-momentum
tensor was suggested by the fact that in the classical 4D $\mathcal{N}=1$ SYM,
the Noether currents associated with the $U(1)_A$ symmetry, SUSY and
translational invariance, form a multiplet under SUSY, the Ferrara--Zumino (FZ)
supermultiplet~\cite{Ferrara:1974pz}.\footnote{See Chapter~20
of~Ref.~\cite{West:1990tg} for a very readable exposition. An interesting
application of the notion of the FZ multiplet in phenomenology was recently
considered in~Ref.~\cite{Kitano:2011fk}.} For the classical
theory~\eqref{eq:(1.4)}, we define the $U(1)_A$
current~$\Breve{\jmath}_{5\mu}(x)$, the SUSY current~$\Breve{S}_\mu(x)$ and the
energy-momentum tensor~$\Breve{T}_{\mu\nu}(x)$ by
\begin{align}
   \Breve{\jmath}_{5\mu}(x)
   &\equiv\tr\left[\Bar{\psi}(x)\gamma_\mu\gamma_5\psi(x)\right],
\notag\\
   \Breve{S}_\mu(x)
   &\equiv-\sigma_{\rho\sigma}\gamma_\mu\tr\left[\psi(x)F_{\rho\sigma}(x)\right],
\notag\\
   \Breve{T}_{\mu\nu}(x)
   &\equiv
   2\tr\left[F_{\mu\rho}(x)F_{\nu\rho}(x)\right]
   -\frac{1}{2}\delta_{\mu\nu}\tr\left[F_{\rho\sigma}(x)F_{\rho\sigma}(x)\right]
\notag\\
   &\qquad{}
   +\frac{1}{4}\tr\left[
   \Bar{\psi}(x)\left(\gamma_\mu\overleftrightarrow{D}_\nu
   +\gamma_\nu\overleftrightarrow{D}_\mu\right)\psi(x)
   \right]
   -\frac{1}{2}\delta_{\mu\nu}\tr\left[
   \Bar{\psi}(x)\overleftrightarrow{\Slash{D}}\psi(x)
   \right],
\label{eq:(1.7)}
\end{align}
where the definitions of~$\Breve{\jmath}_{5\mu}(x)$ and~$\Breve{S}_\mu(x)$ are
standard while, as well-known, the definition of the energy-momentum tensor is
large extent arbitrary. In Eq.~\eqref{eq:(1.7)}, we defined the energy-momentum
tensor by first introducing a background gravitational field
into the continuum action~\eqref{eq:(1.4)} and then taking a flat-space limit
after differentiating the action with respect to the gravitational field. Then
one finds, under the super transformation~\eqref{eq:(1.6)},
\begin{equation}
   \Bar{\delta}_\xi\Breve{\jmath}_{5\mu}(x)=\Bar{\xi}\gamma_5\Breve{S}_\mu(x),
\label{eq:(1.8)}
\end{equation}
and
\begin{align}
   \Bar{\delta}_\xi\Breve{S}_\mu(x)
   &=2\gamma_\nu\xi
   \biggl\{\Breve{T}_{\mu\nu}(x)
   +\frac{3}{4}\delta_{\mu\nu}
   \tr\left[\Bar{\psi}(x)\Slash{D}\psi(x)\right]
\notag\\
   &\qquad\qquad{}
   +\frac{1}{4}\epsilon_{\mu\nu\rho\sigma}
   \partial_\rho\Breve{\jmath}_{5\sigma}(x)
   +\frac{1}{4}\tr\left[\Bar{\psi}(x)\sigma_{\mu\nu}\Slash{D}\psi(x)\right]
   \biggr\}
\notag\\
   &\qquad{}
   +\gamma_5\gamma_\nu\xi\partial_\nu\Breve{\jmath}_{5\mu}(x)
\notag\\
   &\qquad{}
   +\frac{3}{2}\xi\tr\left[\Bar{\psi}(x)\gamma_\mu\Slash{D}\psi(x)\right]
   +\frac{3}{2}\gamma_5\xi\tr\left[\Bar{\psi}(x)\gamma_\mu\gamma_5
   \Slash{D}\psi(x)\right]
\notag\\
   &\qquad{}
   -\frac{1}{2}\gamma_5\gamma_\nu\xi\tr\left[\Bar{\psi}(x)
   \gamma_\mu\gamma_5\gamma_\nu\Slash{D}\psi(x)\right]
\notag\\
   &\qquad{}
   -\frac{1}{4}\sigma_{\nu\rho}\xi\tr\left[\Bar{\psi}(x)
   \gamma_\mu\sigma_{\nu\rho}\Slash{D}\psi(x)\right].
\label{eq:(1.9)}
\end{align}
Similarly, $\Bar{\delta}_\xi\Breve{T}_{\mu\nu}(x)$ becomes a linear combination
of~$\partial_\rho\Breve{S}_\sigma(x)$, up to terms being proportional to the
equation of motion of the (massless) gluino, $\Slash{D}\psi(x)$. By integrating
over the spatial coordinates of the $\mu=0$ component of the above relations,
one sees that they are consistent with the defining algebra between the SUSY
charge, the $U(1)_A$ charge (this is the $R$-charge) and the momentum, up to
the equation of motion $\Slash{D}\psi(x)=0$. Eqs.~\eqref{eq:(1.8)}
and~\eqref{eq:(1.9)} thus may be regarded as a basic characterization of
fundamental Noether currents in~Eq.~\eqref{eq:(1.7)}. By going one step
further, one may regard Eqs.~\eqref{eq:(1.8)} and~\eqref{eq:(1.9)} as defining
relations of the Noether currents. This is the approach we adopt in the present
paper. That is, we define a lattice energy-momentum
tensor~$\mathcal{T}_{\mu\nu}(x)$ by a (renormalized modified) super
transformation of a (renormalized) lattice SUSY current~$\mathcal{S}_\mu(x)$
that is conserved in the continuum limit under an appropriate gluino mass
tuning~\cite{Curci:1986sm}. Then it can be shown that the lattice
energy-momentum tensor so constructed is conserved in the continuum limit. This
conservation of the energy-momentum tensor is not quite trivial, because SUSY
is not a manifest symmetry with the lattice regularization. For the
conservation of~$\mathcal{T}_{\mu\nu}(x)$, a ``finiteness'' of a lattice super
transformation induced by~$\mathcal{S}_\mu(x)$ turns to be crucial. This
finiteness might be regarded as a necessary condition for the existence of a
renormalized 4D $\mathcal{N}=1$ SYM.

As the structure of our energy-momentum tensor is very explicit as
Eqs.~\eqref{eq:(3.3)} and~\eqref{eq:(3.4)} show, it is usable in actual
Monte Carlo simulations of 4D $\mathcal{N}=1$ SYM~\cite{Montvay:1996pz,%
Montvay:1997ak,Koutsoumbas:1997de,Kirchner:1998nk,Kirchner:1998mp,%
Campos:1999du,Feo:1999hw,Feo:1999hx,Farchioni:2000mp,Farchioni:2000kb,%
Farchioni:2001yn,Farchioni:2001yr,Farchioni:2001wx,Peetz:2002sr,%
Farchioni:2004fy,Demmouche:2008ms,Demmouche:2009ki,Demmouche:2010sf,%
Bergner:2011wf,Bergner:2012at,Fleming:2000fa,Giedt:2008xm,Endres:2009yp,%
Endres:2009pu,Kim:2011fw}. Physical quantities such as the viscosity (see, for
example, Ref.~\cite{Sakai:2007cm} and references therein) might be measured
from the two-point function of our lattice energy-momentum tensor.

Our argument below is quite general and rather independent of the details of
a lattice formulation. We assume the locality of the lattice action, that
consists of the sum of a gauge boson (gluon) action~$S_{\text{gluon}}$,
the kinetic term~$S_{\text{gluino}}$ and the mass term~$S_{\text{mass}}$ of the
gauge fermion (gluino),\footnote{Although we use the language of a
four-dimensional lattice theory in the present paper, we believe that our
argument can be transcribed for a
domain-wall-type~\cite{Kaplan:1992bt,Shamir:1993zy} five-dimensional
setting~\cite{Nishimura:1997vg,Maru:1997kh,Kaplan:1999jn} without much
difficulty.}
\begin{equation}
   S\equiv S_{\text{gluon}}+S_{\text{gluino}}+S_{\text{mass}},
\label{eq:(1.10)}
\end{equation}
where
\begin{equation}
   S_{\text{gluino}}\equiv
   a^4\sum_x\tr\left[\Bar{\psi}(x)D\psi(x)\right],
\label{eq:(1.11)}
\end{equation}
where $D$ is a lattice Dirac operator ($(C^{-1}D)^T=-C^{-1}D$), and
\begin{equation}
   S_{\text{mass}}\equiv
   a^4\sum_xM\tr\left[\Bar{\psi}(x)\psi(x)\right].
\label{eq:(1.12)}
\end{equation}
We further assume that the lattice action~\eqref{eq:(1.10)} is invariant under
the hypercubic group~$H(4)$, including the parity transformation~$\mathcal{P}$
defined by
\begin{align}
   U_0(x_0,\Vec{x})&\xrightarrow{\mathcal{P}}U_0(x_0,-\Vec{x}),&
   U_k(x_0,\Vec{x})&\xrightarrow{\mathcal{P}}U_k^\dagger(x_0,-\Vec{x}-a\hat k),
\notag\\
   \psi(x_0,\Vec x)&\xrightarrow{\mathcal{P}}i\gamma_0\psi(x_0,-\Vec x),&
   \Bar{\psi}(x_0,\Vec x)&\xrightarrow{\mathcal{P}}
   -i\Bar{\psi}(x_0,-\Vec x)\gamma_0,
\label{eq:(1.13)}
\end{align}
where and in what follows, $\Vec{x}\equiv(x_1,x_2,x_3)$ and~$k$, $l=1$ or~$2$
or~$3$.

In the following discussion, we make frequent use of a lattice transcription of
the field strength,
\begin{equation}
   \left[F_{\mu\nu}\right]^L(x)\equiv2\tr\left[P_{\mu\nu}(x)T^a\right]T^a,
\label{eq:(1.14)}
\end{equation}
that is defined from the clover plaquette $P_{\mu\nu}(x)$,
\begin{equation}
   P_{\mu\nu}(x)\equiv\frac{1}{4}\sum_{i=1}^4\frac{1}{2ia^2g}
   \left[U_{i\mu\nu}(x)-U_{i\mu\nu}^\dagger(x)\right],
\label{eq:(1.15)}
\end{equation}
where
\begin{align}
   U_{1\mu\nu}(x)
   &\equiv U_\mu(x)U_\nu(x+a\hat\mu)U_\mu^\dagger(x+a\hat\nu)U_\nu^\dagger(x),
\notag\\
   U_{2\mu\nu}(x)
   &\equiv U_\nu(x)U_\mu^\dagger(x-a\hat\mu+a\hat\nu)
   U_\nu^\dagger(x-a\hat\mu)U_\mu(x-a\hat\mu),
\notag\\
   U_{3\mu\nu}(x)
   &\equiv U_\mu^\dagger(x-a\hat\mu)U_\nu^\dagger(x-a\hat\mu-a\hat\nu)
   U_\mu(x-a\hat\mu-a\hat\nu)U_\nu(x-a\hat\nu),
\notag\\
   U_{4\mu\nu}(x)
   &\equiv U_\nu^\dagger(x-a\hat\nu)U_\mu(x-a\hat\nu)
   U_\nu(x+a\hat\mu-a\hat\nu)U_\mu^\dagger(x).
\label{eq:(1.16)}
\end{align}
Note that $[F_{\mu\nu}]^L(x)$ is traceless by construction~\eqref{eq:(1.14)}.
Under the parity, $[F_{\mu\nu}]^L(x)$ transforms in an identical manner as the
continuum field strength. That is, we have
\begin{equation}
   \left[F_{0k}\right]^L(x_0,\Vec{x})
   \xrightarrow{\mathcal{P}}-\left[F_{0k}\right]^L(x_0,-\Vec{x}),\qquad
   \left[F_{kl}\right]^L(x_0,\Vec{x})
   \xrightarrow{\mathcal{P}}+\left[F_{kl}\right]^L(x_0,-\Vec{x}).
\label{eq:(1.17)}
\end{equation}

\section{Lattice SUSY current and a renormalized SUSY WT relation}
\label{sec:2}
In this section, we recall how the proposal of~Ref.~\cite{Curci:1986sm} that
SUSY in the lattice formulation of 4D $\mathcal{N}=1$ SYM is restored with an
appropriate tuning of the gluino mass~$M$ is understood in terms of the SUSY WT
relation~\cite{Curci:1986sm,Taniguchi:1999fc,Farchioni:2001wx,Suzuki:2012pc}.
See also Refs.~\cite{Farchioni:2000kb,Farchioni:2001yr}. This discussion
provides the definition of a lattice SUSY current that is conserved in the
continuum limit and our basic SUSY WT relation that will be fully utilized
in the construction of the energy-momentum tensor in~Sec.~\ref{sec:3}.

First, as a lattice counterpart of the SUSY transformation in the
continuum~\eqref{eq:(1.6)}, we define~\cite{Taniguchi:1999fc}
\begin{align}
   \Bar{\delta}_\xi U_\mu(x)&\equiv iag\frac{1}{2}
   \Bar{\xi}\gamma_\mu\left[
   \psi(x)U_\mu(x)
   +U_\mu(x)\psi(x+a\hat\mu)
   \right],
\notag\\
   \Bar{\delta}_\xi U_\mu^\dagger(x)&\equiv-iag\frac{1}{2}
   \Bar{\xi}\gamma_\mu\left[
   U_\mu^\dagger(x)\psi(x)
   +\psi(x+a\hat\mu)U_\mu^\dagger(x)
   \right],
\notag\\
   \Bar{\delta}_\xi\psi(x)&\equiv-\frac{1}{2}\sigma_{\mu\nu}\xi
   \left[F_{\mu\nu}\right]^L(x),\qquad
   \Bar{\delta}_\xi\Bar{\psi}(x)
   =\frac{1}{2}\Bar{\xi}\sigma_{\mu\nu}\left[F_{\mu\nu}\right]^L(x).
\label{eq:(2.1)}
\end{align}
To derive a WT relation, we need also the localized transformation and,
corresponding to the above, we define
\begin{align}
   \delta_\xi U_\mu(x)&\equiv iag\frac{1}{2}
   \left[
   \Bar{\xi}(x)\gamma_\mu\psi(x)U_\mu(x)
   +\Bar{\xi}(x+a\hat\mu)\gamma_\mu U_\mu(x)\psi(x+a\hat\mu)
   \right],
\notag\\
   \delta_\xi U_\mu^\dagger(x)&\equiv-iag\frac{1}{2}
   \left[
   \Bar{\xi}(x)\gamma_\mu U_\mu^\dagger(x)\psi(x)
   +\Bar{\xi}(x+a\hat\mu)\gamma_\mu\psi(x+a\hat\mu)U_\mu^\dagger(x)
   \right],
\notag\\
   \delta_\xi\psi(x)&\equiv-\frac{1}{2}\sigma_{\mu\nu}\xi(x)
   \left[F_{\mu\nu}\right]^L(x),\qquad
   \delta_\xi\Bar{\psi}(x)
   =\frac{1}{2}\Bar{\xi}(x)\sigma_{\mu\nu}\left[F_{\mu\nu}\right]^L(x),
\label{eq:(2.2)}
\end{align}
where the Grassmann-odd spinor parameter~$\xi(x)$ obeys
$\Bar{\xi}(x)=\xi^T(x)(-C^{-1})$.

Now the lattice action~\eqref{eq:(1.10)} is not invariant
under~$\Bar{\delta}_\xi$~\eqref{eq:(2.1)}, because, first of all, the
gluino mass term~\eqref{eq:(1.12)} explicitly breaks SUSY and, various $O(a)$
lattice artifacts in the lattice action also break SUSY. Since the localized
SUSY transformation~$\delta_\xi$~\eqref{eq:(2.2)} reduces to the global
one~\eqref{eq:(2.1)} for~$\xi(x)\to\xi$, the variation of the lattice
action~\eqref{eq:(1.10)} under~Eq.~\eqref{eq:(2.2)} can be written as
\begin{equation}
   \delta_\xi S=a^4\sum_x\Bar{\xi}(x)
   \left[-\partial_\mu^SS_\mu(x)+M\chi(x)+X_S(x)\right],
\label{eq:(2.3)}
\end{equation}
where combinations $S_\mu(x)$ and~$\chi(x)$ are defined by
\begin{equation}
   S_\mu(x)\equiv-\sigma_{\rho\sigma}\gamma_\mu
   \tr\left\{\psi(x)\left[F_{\rho\sigma}\right]^L(x)\right\},
\label{eq:(2.4)}
\end{equation}
and
\begin{equation}
   \chi(x)\equiv\sigma_{\mu\nu}
   \tr\left\{\psi(x)\left[F_{\mu\nu}\right]^L(x)\right\},
\label{eq:(2.5)}
\end{equation}
respectively. $S_\mu(x)$ in~Eq.~\eqref{eq:(2.4)} is nothing but a lattice
transcription of the SUSY current in the continuum, $\Breve{S}_\mu(x)$
in~Eq.~\eqref{eq:(1.7)}, and $\chi(x)$ in~Eq.~\eqref{eq:(2.5)} is an explicit
SUSY breaking caused by~$S_{\text{mass}}$. Thus, in~Eq.~\eqref{eq:(2.3)},
$X_S(x)$ represents an $O(a)$ SUSY breaking effect attributed to the lattice
artifacts. Thus, considering the infinitesimal variation~\eqref{eq:(2.2)} in
the integration variables of the functional integral, we have an exact identity
that holds for any operator~$\mathcal{O}$,\footnote{Throughout this article, we
always assume that an operator represented by the symbol~$\mathcal{O}$ is gauge
invariant.}
\begin{equation}
   \left\langle\partial_\mu^SS_\mu(x)\mathcal{O}\right\rangle
   =\left\langle\left[M\chi(x)+ X_S(x)\right]\mathcal{O}\right\rangle
   -\left\langle
   \frac{1}{a^4}\frac{\partial}{\partial\Bar{\xi}(x)}\delta_\xi
   \mathcal{O}\right\rangle.
\label{eq:(2.6)}
\end{equation}

Here is a remark on the property of~$X_S(x)$ under the parity: From
Eqs.~\eqref{eq:(2.4)} and~\eqref{eq:(2.5)}, we see that
\begin{equation}
   S_\mu(x_0,\Vec{x})\xrightarrow{\mathcal{P}}
   \begin{cases}
   +i\gamma_0S_0(x_0,-\Vec{x})&\text{for $\mu=0$},\\
   -i\gamma_0S_k(x_0,-\Vec{x})&\text{for $\mu=k$},\\
   \end{cases}
\label{eq:(2.7)}
\end{equation}
and
\begin{equation}
   \chi(x_0,\Vec{x})\xrightarrow{\mathcal{P}}
   i\gamma_0\chi(x_0,-\Vec{x}).
\label{eq:(2.8)}
\end{equation}
These shows that $\partial_\mu^SS_\mu(x)$ possesses the same transformation
property as~$\chi(x)$ under the parity. On the other hand, by defining
$\xi(x_0,\Vec{x})\xrightarrow{\mathcal P}i\gamma_0\xi(x_0,-\Vec{x})$
and~$\Bar{\xi}(x_0,\Vec{x})\xrightarrow{\mathcal P}%
-i\Bar{\xi}(x_0,-\Vec{x})\gamma_0$, we see that $\delta_\xi$
in~Eq.~\eqref{eq:(2.2)} and the parity commute to each other. Then by the
parity invariance of the action~$S$ and Eqs.~\eqref{eq:(2.7)}
and~\eqref{eq:(2.8)}, $X_S(x)$ in~Eq.~\eqref{eq:(2.3)} obeys
\begin{equation}
   X_S(x_0,\Vec{x})\xrightarrow{\mathcal{P}}i\gamma_0X_S(x_0,-\Vec{x}).
\label{eq:(2.9)}
\end{equation}

Next, we consider the renormalization of the composite operator~$X_S(x)$ (just
like for the chiral symmetry~\cite{Bochicchio:1985xa,Testa:1998ez}). It
can be argued that~\cite{Curci:1986sm,Farchioni:2001wx,Suzuki:2012pc}, by using
the assumed symmetries of~$S$ and resultant properties of~$X_S(x)$ (such
as~Eq.~\eqref{eq:(2.9)}), by writing\footnote{In~Ref.~\cite{Suzuki:2012pc},
this operator renormalization was described by using an operator basis in the
continuum theory. With the present lattice regularization, however, it is more
appropriate to use an operator basis in the lattice theory
as~Eq.~\eqref{eq:(2.10)}, because the renormalization is carried out in terms
of counterterms with lattice
structure~\cite{Reisz:1987da,Reisz:1987pw,Luscher:1988sd}. Also, we adopt a
different choice of basis for the term proportional to the equation of motion
from~Ref.~\cite{Suzuki:2012pc}, as this choice is practically more convenient.
Finally, we make use of a simpler lattice SUSY current than
Ref.~\cite{Taniguchi:1999fc}. Then the symmetric difference~$\partial_\mu^S$
instead of the backward difference~$\partial_\mu^*$ appears for the covariance
under the parity.}
\begin{align}
   X_S(x)&=(1-\mathcal{Z}_S)\partial_\mu^SS_\mu(x)
   -\mathcal{Z}_T\partial_\mu^ST_\mu(x)
   -\frac{1}{a}\mathcal{Z}_\chi\chi(x)
\notag\\
   &\qquad{}
   -\mathcal{Z}_{3F}\tr\left[\psi(x)\Bar{\psi}(x)\psi(x)\right]
\notag\\
   &\qquad{}
   -\mathcal{Z}_{\text{EOM}}
   \sigma_{\mu\nu}\tr\{[F_{\mu\nu}]^L(x)(D+M)\psi(x)\}
\notag\\
   &\qquad{}
   +a\mathcal{E}(x),
\label{eq:(2.10)}
\end{align}
and choosing renormalization constants $\mathcal{Z}_S$, $\mathcal{Z}_T$,
$\mathcal{Z}_\chi$, $\mathcal{Z}_{3F}$ and~$\mathcal{Z}_{\text{EOM}}$
appropriately, the dimension~$11/2$ operator $\mathcal{E}(x)$ can be made at
most logarithmically divergent for~$a\to0$. More precisely, the operator
$\mathcal{E}(x)$ is given by a linear combination of renormalized operators
with logarithmically divergent coefficients. Since the operator~$\mathcal{O}$
in~Eq.~\eqref{eq:(2.6)} is gauge invariant, the mixing of~$X_S(x)$ with gauge
non-invariant operators~\cite{Taniguchi:1999fc} does not occur
in~Eq.~\eqref{eq:(2.10)}. A new lattice operator~$T_\mu(x)$
in~Eq.~\eqref{eq:(2.10)} is defined by
\begin{equation}
   T_\mu(x)\equiv2\gamma_\nu
   \tr\left\{\psi(x)\left[F_{\mu\nu}\right]^L(x)\right\},
\label{eq:(2.11)}
\end{equation}
and, from this definition,
\begin{equation}
   T_\mu(x_0,\Vec{x})\xrightarrow{\mathcal{P}}
   \begin{cases}
   +i\gamma_0T_0(x_0,-\Vec{x})&\text{for $\mu=0$},\\
   -i\gamma_0T_k(x_0,-\Vec{x})&\text{for $\mu=k$}.\\
   \end{cases}
\label{eq:(2.12)}
\end{equation}

From Eqs.~\eqref{eq:(2.6)} and~\eqref{eq:(2.10)}, we have
\begin{align}
   &\left\langle
   \partial_\mu^S
   \left[\mathcal{Z}_SS_\mu(x)+\mathcal{Z}_TT_\mu(x)\right]
   \mathcal{O}\right\rangle
\notag\\
   &=\left(M-\frac{1}{a}\mathcal{Z}_\chi\right)
   \left\langle
   \chi(x)\mathcal{O}\right\rangle
   -\mathcal{Z}_{3F}
   \left\langle\tr\left[\psi(x)\Bar{\psi}(x)\psi(x)\right]
   \mathcal{O}\right\rangle
\notag\\
   &\qquad{}
   -\left\langle
   \frac{1}{a^4}\frac{\partial}{\partial\Bar{\xi}(x)}\delta_\xi
   \mathcal{O}\right\rangle
   -\mathcal{Z}_{\text{EOM}}\left\langle
   \sigma_{\mu\nu}\tr\{[F_{\mu\nu}]^L(x)(D+M)\psi(x)\}
   \mathcal{O}\right\rangle
\notag\\
   &\qquad\qquad{}+\left\langle a\mathcal{E}(x)
   \mathcal{O}\right\rangle.
\label{eq:(2.13)}
\end{align}
In this expression, we first note
\begin{equation}
   \left\langle
   \sigma_{\mu\nu}\tr\{[F_{\mu\nu}]^L(x)(D+M)\psi(x)\}
   \mathcal{O}\right\rangle
   =\left\langle
   \frac{1}{a^4}\frac{\partial}{\partial\Bar{\xi}(x)}\delta_{F\xi}
   \mathcal{O}\right\rangle,
\label{eq:(2.14)}
\end{equation}
where the localized transformation~$\delta_{F\xi}$ is a super-like
transformation that acts non-trivially only on the gluino:
\begin{equation}
   \delta_{F\xi}U_\mu(x)=0,\qquad
   \delta_{F\xi}\psi(x)=\delta_\xi\psi(x),\qquad
   \delta_{F\xi}\Bar{\psi}(x)=\delta_\xi\Bar{\psi}(x),
\label{eq:(2.15)}
\end{equation}
because the left-hand side of~Eq.~\eqref{eq:(2.14)} is proportional to the
equation of motion of the gluino (i.e., this is the Schwinger--Dyson (SD)
equation). Therefore, we have
\begin{align}
   &\left\langle
   \partial_\mu^S
   \left[\mathcal{Z}_SS_\mu(x)+\mathcal{Z}_TT_\mu(x)\right]
   \mathcal{O}\right\rangle
\notag\\
   &\qquad=\left(M-\frac{1}{a}\mathcal{Z}_\chi\right)
   \left\langle
   \chi(x)\mathcal{O}\right\rangle
   -\mathcal{Z}_{3F}
   \left\langle\tr\left[\psi(x)\Bar{\psi}(x)\psi(x)\right]
   \mathcal{O}\right\rangle
\notag\\
   &\qquad\qquad{}
   +\left\langle
   \left[-\frac{1}{a^4}\frac{\partial}{\partial\Bar{\xi}(x)}\Delta_\xi
   +a\mathcal{E}(x)\right]
   \mathcal{O}\right\rangle,
\label{eq:(2.16)}
\end{align}
where $\Delta_\xi$ is a \emph{modified\/} localized super transformation,
\begin{equation}
   \Delta_\xi\equiv
   \delta_\xi+\mathcal{Z}_{\text{EOM}}\delta_{F\xi}.
\label{eq:(2.17)}
\end{equation}

Let us consider a special case of~Eq.~\eqref{eq:(2.16)} that $\mathcal{O}$
is a collection of renormalized local operators and the point~$x$ stays away
from the support of~$\mathcal{O}$ by a finite physical distance. In what
follows, we express this situation by
\begin{equation}
   x\leftrightsquigarrow\supp(\mathcal{O}).
\label{eq:(2.18)}
\end{equation}
Under this situation, the $\Bar{\xi}(x)$-derivative in the last line
of~Eq.~\eqref{eq:(2.16)} identically vanishes.\footnote{Our argument below
works when the operator~$\mathcal{O}$ is ultra-local, i.e., its support is
a strictly-finite region on the lattice. With the use of the overlap lattice
Dirac operator~$D$~\cite{Neuberger:1997fp,Neuberger:1998wv}, however,
one might be interested in an operator~$\mathcal{O}$ that is not ultra-local
but exponentially local~\cite{Hernandez:1998et}. Although we do not give any
analysis for such a case, on physical grounds, we believe that our conclusions
will not change.} Also the last term of~Eq.~\eqref{eq:(2.16)} vanishes in the
continuum limit, because for~$x\leftrightsquigarrow\supp(\mathcal{O})$ the
dimension~$11/2$ operator~$\mathcal{E}(x)$ does not produce any $O(1/a)$
divergence that can compensate the overall factor of~$a$. We thus have
\begin{align}
   &\left\langle
   \partial_\mu^S
   \left[\mathcal{Z}_SS_\mu(x)+\mathcal{Z}_TT_\mu(x)\right]
   \mathcal{O}\right\rangle
\notag\\
   &\qquad{}\xrightarrow{a\to0}
   \left(M-\frac{1}{a}\mathcal{Z}_\chi\right)
   \left\langle
   \chi(x)\mathcal{O}\right\rangle
   -\mathcal{Z}_{3F}
   \left\langle\tr\left[\psi(x)\Bar{\psi}(x)\psi(x)\right]
   \mathcal{O}\right\rangle,
\notag\\
   &\qquad\qquad\text{for $x\leftrightsquigarrow\supp(\mathcal{O})$}.
\label{eq:(2.19)}
\end{align}
This can be regarded as a would-be conservation law of a lattice SUSY current.
If this relation with~$\mathcal{Z}_{\text{3F}}\neq0$ held, then the last term,
that is cubic in the gluino field, would give rise to an ``exotic'' SUSY
anomaly. This is not what we expect and we thus assume its absence:
\begin{equation}
   \mathcal{Z}_{3F}=0.
\label{eq:(2.20)}
\end{equation}
As shown in~Ref.~\cite{Suzuki:2012pc} by utilizing the generalized BRS
transformation, this is actually the case at least to all orders in the
perturbation theory. Accepting Eq.~\eqref{eq:(2.20)}, we then have
\begin{align}
   &\left\langle
   \partial_\mu^S
   \left[\mathcal{Z}_SS_\mu(x)+\mathcal{Z}_TT_\mu(x)\right]
   \mathcal{O}\right\rangle
   \xrightarrow{a\to0}
   \left(M-\frac{1}{a}\mathcal{Z}_\chi\right)
   \left\langle
   \chi(x)\mathcal{O}\right\rangle,
\notag\\
   &\qquad{}\text{for $x\leftrightsquigarrow\supp(\mathcal{O})$}.
\label{eq:(2.21)}
\end{align}
The combination $M-(1/a)\mathcal{Z}_\chi$ in the right-hand side corresponds to
an additive renormalization of the gluino mass. So we tune the bare mass
parameter~$M$~\cite{Curci:1986sm} so that
\begin{equation}
   M-\frac{1}{a}\mathcal{Z}_\chi=0.
\label{eq:(2.22)}
\end{equation}
In actual numerical simulations, this tuning should be carried out for each
value of the lattice spacing~$a$. In this way, we obtain the conservation
law in the continuum limit:
\begin{equation}
   \left\langle
   \partial_\mu^S
   \left[\mathcal{Z}_SS_\mu(x)+\mathcal{Z}_TT_\mu(x)\right]
   \mathcal{O}\right\rangle
   \xrightarrow{a\to0}0,
   \qquad\text{for $x\leftrightsquigarrow\supp(\mathcal{O})$},
\label{eq:(2.23)}
\end{equation}
that is a minimal requirement for the restoration of SUSY in the continuum
limit.

We consider the renormalization of the conserved current
$\mathcal{Z}_SS_\mu(x)+\mathcal{Z}_TT_\mu(x)$. Since this is a dimension~$7/2$
gauge invariant vectorial operator, it can mix only with the operators
$S_\mu(x)$ and~$T_\mu(x)$. However, if the mixing were generic, it would be
inconsistent with the conservation law~\eqref{eq:(2.23)}. That is, the
renormalization must be multiplicative.\footnote{It is very interesting to
confirm this renormalizability by the lattice perturbation theory along the
line of~Ref.~\cite{Taniguchi:1999fc}.} Thus we set
\begin{equation}
   \mathcal{S}_\mu(x)
   \equiv\mathcal{Z}
   \left[\mathcal{Z}_SS_\mu(x)+\mathcal{Z}_TT_\mu(x)\right],
\label{eq:(2.24)}
\end{equation}
and call this a renormalized lattice SUSY current. With an appropriate
choice of the constant~$\mathcal{Z}$, $\mathcal{S}_\mu(x)$ should have finite
correlation functions with any renormalized operator~$\mathcal{O}$, as
far as~$x\leftrightsquigarrow\supp(\mathcal{O})$. The renormalization constants
$\mathcal{Z}_S$ and~$\mathcal{Z}_T$ in~Eqs.~\eqref{eq:(2.10)}
and~\eqref{eq:(2.24)} emerge from the power-divergence subtraction and it can
be argued that such constants are independent of the renormalization
scale~$\mu^2$ and finite as~$a\to0$~\cite{Testa:1998ez}. On the other hand,
by the dimensional reason the constant~$\mathcal{Z}$ in~Eq.~\eqref{eq:(2.24)}
is at most logarithmically divergent. In terms of the renormalized SUSY
current~\eqref{eq:(2.24)}, the SUSY WT relation~\eqref{eq:(2.16)} under the
assumptions~\eqref{eq:(2.20)} and~\eqref{eq:(2.22)} reads
\begin{equation}
   \left\langle
   \partial_\mu^S\mathcal{S}_\mu(x)
   \mathcal{O}\right\rangle
   =\left\langle
   \mathcal{Z}
   \left[-\frac{1}{a^4}\frac{\partial}{\partial\Bar{\xi}(x)}\Delta_\xi
   +a\mathcal{E}(x)\right]
   \mathcal{O}\right\rangle.
\label{eq:(2.25)}
\end{equation}

Finally, we argue that for any renormalized local operator~$\mathcal{O}$, the
combination
\begin{equation}
   \mathcal{Z}
   \left[-\frac{1}{a^4}\frac{\partial}{\partial\Bar{\xi}(x)}\Delta_\xi
   +a\mathcal{E}(x)\right]
   \mathcal{O},
\label{eq:(2.26)}
\end{equation}
when $x\in\supp(\mathcal{O})$, produces a renormalized local operator in the
support of~$\mathcal{O}$. In~Eq.~\eqref{eq:(2.26)}, something non-trivial in
the continuum limit can occur only when $x\in\supp(\mathcal{O})$, because
otherwise the $\Bar{\xi}(x)$-derivative vanishes and the dimension~$11/2$
operator~$\mathcal{E}(x)$ does not produce any $O(1/a)$ divergence that can
compensate the overall factor of~$a$. To see what happens when
$x\in\supp(\mathcal{O})$, we sum Eq.~\eqref{eq:(2.25)} over the point~$x$
within a finite region~$\mathcal{D}_{\mathcal{O}}$ that contains the
operator~$\mathcal{O}$. Then, since
$\sum_{x\in\mathcal{D}_{\mathcal{O}}}\partial_\mu^S\mathcal{S}_\mu(x)$ in the
left-hand side becomes a collection of (mutually non-overlapping) renormalized
operators that have no overlap with~$\mathcal{O}$, the left-hand side remains
finite as $a\to0$. This shows that the sum~$\sum_{x\in\mathcal{D}_{\mathcal{O}}}$
of the right-hand side of~Eq.~\eqref{eq:(2.25)} is finite too and we see that
the combination~\eqref{eq:(2.26)} does not produce any ultraviolet divergence
even for $x\in\supp(\mathcal{O})$. We note that such a finiteness will be
indispensable for the existence of a renormalized SYM, because for that there
should exist a renormalized SUSY current that generates a finite super
transformation on renormalized fields. Since the above analysis shows that
$\mathcal{S}_\mu(x)$ is the unique conserved SUSY current in the present
framework, $\mathcal{S}_\mu(x)$ should generate such a finite super
transformation. This argument on the basis of the existence of a renormalized
SYM also supports the finiteness of the combination~\eqref{eq:(2.26)}.

\section{A lattice energy-momentum tensor and its conservation law}
\label{sec:3}
We now define an energy-momentum tensor~$\mathcal{T}_{\mu\nu}(x)$ on the lattice
through the relation quite analogous to~Eq.~\eqref{eq:(1.9)}. That is,
\begin{align}
   \mathcal{Z}\Bar{\Delta}_\xi\mathcal{S}_\mu(x)
   &\equiv
   2\gamma_\nu\xi\bigl\{
   \mathcal{T}_{\mu\nu}(x)+c\delta_{\mu\nu}
   \tr\left[\Bar{\psi}(x)(D+M)\psi(x)\right]
\notag\\
   &\qquad\qquad{}
   +(\text{terms anti-symmetric in $\mu$ and~$\nu$})\bigr\}
\notag\\
   &\qquad{}
   +(\text{terms proportional to $\gamma_5\gamma_\nu\xi$, $\xi$,
   $\gamma_5\xi$, $\sigma_{\nu\rho}\xi$}),
\label{eq:(3.1)}
\end{align}
where $c$~is a constant and $\Bar{\Delta}_\xi$ is a ``global version'' of the
localized super transformation~\eqref{eq:(2.17)}, i.e.,
\begin{equation}
   \Bar{\Delta}_\xi\equiv\Bar{\delta}_\xi
   +\mathcal{Z}_{\text{EOM}}\Bar{\delta}_{F\xi},
\label{eq:(3.2)}
\end{equation}
and $\Bar{\delta}_{F\xi}$ is obtained by setting~$\xi(x)\to\xi$
in~Eq.~\eqref{eq:(2.15)}. One can invert Eq.~\eqref{eq:(3.1)} with respect
to~$\mathcal{T}_{\mu\nu}(x)$: We introduce (here and in what follows, indices
$\alpha$, $\beta$, \dots, which run over $1$, $2$, $3$, $4$, denote the spinor
indices),
\begin{align}
   \varTheta_{\mu\nu}(x)
   &\equiv\frac{1}{8}(\gamma_\nu)_{\beta\alpha}
   \frac{\partial}{\partial\xi_\beta}
   \left[\mathcal{Z}\Bar{\Delta}_\xi\mathcal{S}_\mu(x)\right]_\alpha
\notag\\
   &=\mathcal{Z}^2\mathcal{Z}_S
   \frac{1}{8}(\gamma_\nu)_{\beta\alpha}
   \frac{\partial}{\partial\xi_\beta}
   \left\{\left(\Bar{\delta}_\xi
   +\mathcal{Z}_{\text{EOM}}\Bar{\delta}_{F\xi}\right)\left[
   S_\mu(x)+\frac{\mathcal{Z}_T}{\mathcal{Z}_S}T_\mu(x)
   \right]\right\}_\alpha,
\label{eq:(3.3)}
\end{align}
where in the second equality we have substituted Eqs.~\eqref{eq:(3.2)}
and~\eqref{eq:(2.24)}, then
\begin{equation}
   \mathcal{T}_{\mu\nu}(x)
   =\frac{1}{2}\left[\varTheta_{\mu\nu}(x)+\varTheta_{\nu\mu}(x)\right]
   -c\delta_{\mu\nu}
   \tr\left[\Bar{\psi}(x)(D+M)\psi(x)\right].
\label{eq:(3.4)}
\end{equation}

From relations obtained so far, one can confirm that $\varTheta_{\mu\nu}(x)$
in~Eq.~\eqref{eq:(3.3)} behaves as a second-rank tensor under the
parity:\footnote{To see this, it is convenient to note that $\Bar{\Delta}_\xi$
and the parity commute to each other, if one assigns
$\xi\xrightarrow{\mathcal{P}}i\gamma_0\xi$ and
$\Bar{\xi}\xrightarrow{\mathcal{P}}-i\Bar{\xi}\gamma_0$.}
\begin{align}
   \varTheta_{\mu\nu}(x_0,\Vec{x})&\xrightarrow{\mathcal{P}}
   \begin{cases}
   +\varTheta_{00}(x_0,-\Vec{x})&\text{for $\mu=\nu=0$},\\
   -\varTheta_{k0}(x_0,-\Vec{x})&\text{for $\mu=k$ and $\nu=0$},\\
   -\varTheta_{0k}(x_0,-\Vec{x})&\text{for $\mu=0$ and $\nu=k$},\\
   +\varTheta_{kl}(x_0,-\Vec{x})&\text{for $\mu=k$ and $\nu=l$}.
   \end{cases}
\label{eq:(3.5)}
\end{align}

The defining relation~\eqref{eq:(3.1)} is quite analogous to the classical
relation~\eqref{eq:(1.9)}. Yet they differ in that the super
transformation~$\Bar{\delta}_\xi$ is replaced by the renormalized modified
super transformation~$\mathcal{Z}\Bar{\Delta}_\xi$ and the Dirac operator is
shifted by the gluino mass~$M$; also the coefficient of the term that is
proportional to the fermion action is changed from~$3/4$ to~$c$.

Quite interestingly, we can show that the energy-momentum tensor defined
by~Eq.~\eqref{eq:(3.4)} is conserved in the quantum continuum limit. That is
\begin{equation}
   \left\langle\partial_\mu^S\mathcal{T}_{\mu\nu}(x)\mathcal{O}\right\rangle
   \xrightarrow{a\to0}0,
   \qquad\text{for $x\leftrightsquigarrow\supp(\mathcal{O})$},
\label{eq:(3.6)}
\end{equation}
for any renormalized local operator~$\mathcal{O}$. Note that this conservation
law holds for any value of~$c$ in~Eq.~\eqref{eq:(3.4)} because the term
$\tr[\Bar{\psi}(x)(D+M)\psi(x)]$ is proportional to the equation of motion and
it has no correlation with~$\mathcal{O}$ when $x$ is not in the support
of~$\mathcal{O}$. When $x\notin\supp(\mathcal{O})$, one can replace it by
\begin{equation}
   \left\langle\tr\left[\Bar{\psi}(x)(D+M)\psi(x)\right]\right\rangle
   =-2(N_c^2-1)a^{-4}.
\label{eq:(3.7)}
\end{equation}
Since this is independent of~$x$, the last term in~Eq.~\eqref{eq:(3.4)} does not
affect the conservation law. The most general structure of the energy-momentum
tensor (that is symmetric in its indices) in the lattice formulation of the
Yang--Mills theory coupled to fermions has been given
in~Ref.~\cite{Caracciolo:1989pt}; see also
Refs.~\cite{Caracciolo:1989bu,Caracciolo:1991vc,Caracciolo:1991cp} for
related analyses. According to~Ref.~\cite{Caracciolo:1989pt}, a symmetric
energy-momentum tensor that satisfies the conservation law in the continuum
limit~\eqref{eq:(3.6)}, if it exists, is essentially unique. Only ambiguities
are the overall normalization ($\mathcal{Z}^2\mathcal{Z}_S$
in~Eq.~\eqref{eq:(3.3)} in our construction) and a proportionality constant to
the fermion Lagrangian (i.e., the constant $c$ in~Eq.~\eqref{eq:(3.4)}). The
bottom line is that our definition~\eqref{eq:(3.4)}, that is suggested from the
structure of the FZ supermultiplet, gives rise to a very explicit form of this
unique conserved symmetric energy-momentum tensor.

To show the conservation law~\eqref{eq:(3.6)}, we first show the conservation
law of~$\varTheta_{\mu\nu}(x)$ in~Eq.~\eqref{eq:(3.3)} that is not necessary
symmetric under~$\mu\leftrightarrow\nu$. Then we show the anti-symmetric
part of~$\varTheta_{\mu\nu}(x)$ also is conserved in the continuum limit,
implying~Eq.~\eqref{eq:(3.6)}.

To show the conservation law of~$\varTheta_{\mu\nu}(x)$, we set
$\mathcal{O}\to\partial_\nu^S\mathcal{S}_\nu(y)\mathcal{O}$ in the WT
relation~\eqref{eq:(2.25)}. This yields
\begin{align}
   &\left\langle
   \left[\partial_\mu^S\mathcal{S}_\mu(x)\right]_\alpha
   \left[\partial_\nu^S\mathcal{S}_\nu(y)\right]_\beta
   \mathcal{O}\right\rangle
\notag\\
   &=\left\langle
   \mathcal{Z}\left[
   -\frac{1}{a^4}\frac{\partial}{\partial\Bar{\xi}(x)}\Delta_\xi
   +a\mathcal{E}(x)\right]_\alpha
   \left[\partial_\nu^S\mathcal{S}_\nu(y)\right]_\beta
   \mathcal{O}\right\rangle
\notag\\
   &=-\left\langle
   \mathcal{Z}
   \left[
   -\frac{1}{a^4}\frac{\partial}{\partial\Bar{\xi}(y)}\Delta_\xi
   +a\mathcal{E}(y)\right]_\beta
   \left[\partial_\mu^S\mathcal{S}_\mu(x)\right]_\alpha
   \mathcal{O}\right\rangle
\notag\\
   &=\left\langle
   \frac{1}{a^4}\frac{\partial}{\partial\Bar{\xi}_\beta(y)}
   \left[\mathcal{Z}\Delta_\xi\partial_\mu^S\mathcal{S}_\mu(x)\right]_\alpha
   \mathcal{O}\right\rangle
\notag\\
   &\qquad{}+\left\langle
   \left[\partial_\mu^S\mathcal{S}_\mu(x)\right]_\alpha
   \mathcal{Z}\left[
   -\frac{1}{a^4}\frac{\partial}{\partial\Bar{\xi}(y)}\Delta_\xi
   +a\mathcal{E}(y)\right]_\beta
   \mathcal{O}\right\rangle.
\label{eq:(3.8)}
\end{align}
The second equality holds because the first expression is anti-symmetric under
the exchange, $x\leftrightarrow y$, $\mu\leftrightarrow\nu$
and~$\alpha\leftrightarrow\beta$. Using Eq.~\eqref{eq:(2.25)} once again in
the last line, we have the identity
\begin{align}
   &\left\langle
   \frac{1}{a^4}\frac{\partial}{\partial\Bar{\xi}_\beta(y)}
   \left[\mathcal{Z}\Delta_\xi\partial_\mu^S\mathcal{S}_\mu(x)\right]_\alpha
   \mathcal{O}\right\rangle
\notag\\
   &=\left\langle
   \mathcal{Z}\left[
   -\frac{1}{a^4}\frac{\partial}{\partial\Bar{\xi}(x)}
   \Delta_\xi+a\mathcal{E}(x)\right]_\alpha
   \left[\partial_\nu^S\mathcal{S}_\nu(y)\right]_\beta
   \mathcal{O}\right\rangle
\notag\\
   &\qquad{}-\left\langle
   \mathcal{Z}\left[
   -\frac{1}{a^4}\frac{\partial}{\partial\Bar{\xi}(x)}
   \Delta_\xi+a\mathcal{E}(x)\right]_\alpha
   \mathcal{Z}\left[
   -\frac{1}{a^4}\frac{\partial}{\partial\Bar{\xi}(y)}
   \Delta_\xi+a\mathcal{E}(y)\right]_\beta
   \mathcal{O}\right\rangle.
\label{eq:(3.9)}
\end{align}
We then sum this relation over~$y$ within a finite region~$\mathcal{D}_x$
containing the operator~$\partial_\mu^S\mathcal{S}_\mu(x)$. Noting the
identity between the local and global transformations,
\begin{equation}
   \sum_{y\in\mathcal{D}_x}\frac{\partial}{\partial\Bar{\xi}_\beta(y)}\Delta_\xi
   \partial_\mu^S\mathcal{S}_\mu(x)
   =\frac{\partial}{\partial\Bar{\xi}_\beta}\Bar{\Delta}_\xi
   \partial_\mu^S\mathcal{S}_\mu(x),
\label{eq:(3.10))}
\end{equation}
we have
\begin{align}
   &\left\langle
   \frac{\partial}{\partial\Bar{\xi}_\beta}
   \left[\mathcal{Z}\Bar{\Delta}_\xi
   \partial_\mu^S\mathcal{S}_\mu(x)\right]_\alpha
   \mathcal{O}\right\rangle
\notag\\
   &=\left\langle
   \mathcal{Z}\left[
   -\frac{1}{a^4}\frac{\partial}{\partial\Bar{\xi}(x)}
   \Delta_\xi+a\mathcal{E}(x)\right]_\alpha
   a^4\sum_{y\in\mathcal{D}_x}\left[\partial_\nu^S\mathcal{S}_\nu(y)\right]_\beta
   \mathcal{O}\right\rangle
\notag\\
   &\qquad{}-\left\langle
   \mathcal{Z}\left[
   -\frac{1}{a^4}\frac{\partial}{\partial\Bar{\xi}(x)}\Delta_\xi
   +a\mathcal{E}(x)\right]_\alpha
   a^4\sum_{y\in\mathcal{D}_x}\mathcal{Z}\left[
   -\frac{1}{a^4}\frac{\partial}{\partial\Bar{\xi}(y)}\Delta_\xi
   +a\mathcal{E}(y)\right]_\beta
   \mathcal{O}\right\rangle.
\label{eq:(3.11)}
\end{align}
Since $\Bar{\xi}=\xi^T(-C^{-1})$ implies
\begin{equation}
   (\gamma_\nu)_{\beta\alpha}
   \frac{\partial}{\partial\xi_\beta}
   \left[\mathcal{Z}\Bar{\Delta}_\xi
   \partial_\mu^S\mathcal{S}_\mu(x)\right]_\alpha
   =(C^{-1}\gamma_\nu)_{\alpha\beta}
   \frac{\partial}{\partial\Bar{\xi}_\beta}
   \left[\mathcal{Z}\Bar{\Delta}_\xi
   \partial_\mu^S\mathcal{S}_\mu(x)\right]_\alpha,
\label{eq:(3.12)}
\end{equation}
Eqs.~\eqref{eq:(3.3)} and~\eqref{eq:(3.11)} yield,
\begin{align}
   &\left\langle
   \partial_\mu^S\varTheta_{\mu\nu}(x)
   \mathcal{O}\right\rangle
\notag\\
   &=\frac{1}{8}
   (C^{-1}\gamma_\nu)_{\alpha\beta}
   \Biggl\langle
   \mathcal{Z}\left[
   -\frac{1}{a^4}\frac{\partial}{\partial\Bar{\xi}(x)}\Delta_\xi
   +a\mathcal{E}(x)\right]_\alpha
   a^4\sum_{y\in\mathcal{D}_x}
   \left[\partial_\nu^S\mathcal{S}_\nu(y)\right]_\beta
   \mathcal{O}\Biggr\rangle
\notag\\
   &\qquad{}
   -\frac{1}{8}
   (C^{-1}\gamma_\nu)_{\alpha\beta}
   \Biggl\langle
   \mathcal{Z}\left[
   -\frac{1}{a^4}\frac{\partial}{\partial\Bar{\xi}(x)}\Delta_\xi
   +a\mathcal{E}(x)\right]_\alpha
\notag\\
   &\qquad\qquad\qquad\qquad\qquad{}
   \times a^4\sum_{y\in\mathcal{D}_x}\mathcal{Z}\left[
   -\frac{1}{a^4}\frac{\partial}{\partial\Bar{\xi}(y)}
   \Delta_\xi+a\mathcal{E}(y)\right]_\beta
   \mathcal{O}\Biggr\rangle.
\label{eq:(3.13)}
\end{align}
Now, let us suppose that the point~$x$ stays away from the support
of the renormalized operator~$\mathcal{O}$ by a finite physical distance, i.e.,
$x\leftrightsquigarrow\supp(\mathcal{O})$. Let us further assume that the
region~$\mathcal{D}_x$ is chosen so that it does not overlap with the support
of~$\mathcal{O}$. In this situation, since
$\sum_{y\in\mathcal{D}_x}\partial_\nu^S\mathcal{S}_\nu(y)$ does not have any
support at~$x$, Eq.~\eqref{eq:(3.13)} reduces to
\begin{align}
   &\left\langle
   \partial_\mu^S\varTheta_{\mu\nu}(x)
   \mathcal{O}\right\rangle
\notag\\
   &=\frac{1}{8}
   (C^{-1}\gamma_\nu)_{\alpha\beta}
   \Biggl\langle
   \mathcal{Z}\left[
   a\mathcal{E}(x)\right]_\alpha
   a^4\sum_{y\in\mathcal{D}_x}
   \left[\partial_\nu^S\mathcal{S}_\nu(y)\right]_\beta
   \mathcal{O}\Biggr\rangle
\notag\\
   &\qquad{}
   -\frac{1}{8}
   (C^{-1}\gamma_\nu)_{\alpha\beta}
   \Biggl\langle
    a^4\sum_{y\in\mathcal{D}_x}
   \mathcal{Z}\left[
   -\frac{1}{a^4}\frac{\partial}{\partial\Bar{\xi}(x)}\Delta_\xi
   +a\mathcal{E}(x)\right]_\alpha
   \mathcal{Z}\left[a\mathcal{E}(y)\right]_\beta
   \mathcal{O}\Biggr\rangle.
\label{eq:(3.14)}
\end{align}
The first term in the right-hand side is a correlation function of renormalized
operators without any mutual overlap with the overall factor of~$a$. Thus, this
term vanishes in the $a\to0$ limit. In the second term in the right-hand side,
we use the property argued around Eq.~\eqref{eq:(2.26)}. Then it is also a
correlation function of renormalized operators without any mutual overlap with
the overall factor of~$a$. Thus this also vanishes in the $a\to0$ limit. In
this way, we conclude
\begin{equation}
   \left\langle\partial_\mu^S\varTheta_{\mu\nu}(x)\mathcal{O}\right\rangle
   \xrightarrow{a\to0}0,
   \qquad\text{for $x\leftrightsquigarrow\supp(\mathcal{O})$},
\label{eq:(3.15)}
\end{equation}
for any renormalized operator~$\mathcal{O}$.

Next, we consider the anti-symmetric part of~$\varTheta_{\mu\nu}(x)$,
\begin{equation}
   \mathcal{A}_{\mu\nu}(x)
   \equiv\frac{1}{2}
   \left[\varTheta_{\mu\nu}(x)-\varTheta_{\nu\mu}(x)\right],
\label{eq:(3.16)}
\end{equation}
and its renormalization. It must be possible to expand this composite operator
by gauge invariant operators with dimensions less than or equal to~$4$ that
behave in the same way under the hypercubic group as the right-hand side
of~Eq.~\eqref{eq:(3.16)} (for the parity, see~Eq.~\eqref{eq:(3.5)}). By taking
also the constraint~\eqref{eq:(1.5)} into account, the most general possibility
turns to be\footnote{For example, a seemingly-obvious candidate
$\tr[\Bar{\psi}(x)\sigma_{\mu\nu}\psi(x)]$ identically vanishes because
of~Eq.~\eqref{eq:(1.5)}; another candidate
$\tr[\Bar{\psi}(x)\gamma_\mu D_\nu^S\psi(x)]-(\mu\leftrightarrow\nu)$, where
$D_\mu^Sf(x)\equiv(1/2a)[U_\mu(x)f(x+a\Hat{\mu})U_\mu^\dagger(x)
-U_\mu(x-a\Hat{\mu})^\dagger f(x-a\Hat{\mu})U_\mu(x-a\Hat{\mu})]$,
reduces to~Eq.~\eqref{eq:(3.17)}, by the fact that $D=\Slash{D}+O(a)$
and~$\sigma_{\mu\nu}\Slash{D}
=-\epsilon_{\mu\nu\rho\sigma}D_\rho\gamma_\sigma\gamma_5
+\gamma_\mu D_\nu-\gamma_\nu D_\mu$.}
\begin{align}
   &\mathcal{A}_{\mu\nu}(x)
\notag\\
   &=A_1\epsilon_{\mu\nu\rho\sigma}\partial_\rho^S
   \tr\left[\Bar{\psi}(x)\gamma_\sigma\gamma_5\psi(x)\right]
   +A_2\tr\left[\Bar{\psi}(x)\sigma_{\mu\nu}(D+M)\psi(x)\right]
   +a\mathcal{G}_{\mu\nu}(x).
\label{eq:(3.17)}
\end{align}
That is, by choosing the constants $A_1$ and~$A_2$ appropriately, the
dimension~$5$ operator~$\mathcal{G}_{\mu\nu}(x)$ can be made at most
logarithmically divergent. Now, quite fortunately, Eq.~\eqref{eq:(3.17)} shows
that the anti-symmetric part is conserved by itself:
\begin{equation}
   \left\langle\partial_\mu^S\mathcal{A}_{\mu\nu}(x)\mathcal{O}\right\rangle
   \xrightarrow{a\to0}0,
   \qquad\text{for $x\leftrightsquigarrow\supp(\mathcal{O})$}.
\label{eq:(3.18)}
\end{equation}
This is trivially true for the first term of~Eq.~\eqref{eq:(3.17)}. For the
second term, this follows from the SD equation (the equation of motion
of~$\psi(x)$). Finally, the last term of~Eq.~\eqref{eq:(3.17)} vanishes in
the continuum limit. Combined Eqs.~\eqref{eq:(3.18)} and~\eqref{eq:(3.15)},
we conclude that the symmetric part of~$\varTheta_{\mu\nu}(x)$,
$(1/2)[\varTheta_{\mu\nu}(x)+\varTheta_{\nu\mu}(x)]$, also is conserved. This
proves our assertion Eq.~\eqref{eq:(3.6)}, because the symmetric part
of~$\varTheta_{\mu\nu}(x)$ is the first term of~Eq.~\eqref{eq:(3.4)} and the
term with the coefficient~$c$ does not affect the conservation law as we
already explained.

We note that, with our definition of the energy-momentum
tensor~\eqref{eq:(3.4)}, the trace
anomaly~\cite{Crewther:1972kn,Chanowitz:1972vd,Adler:1976zt,Nielsen:1977sy}
and the gamma-trace anomaly (superconformal
anomaly)~\cite{Abbott:1977in,Curtright:1977cg,Inagaki:1978iu,Majumdar:1980ej,
Nicolai:1980km,Hagiwara:1979pu,Hagiwara:1980ys,Kumar:1982ng,Nakayama:1983qt}
are related by\footnote{Ref.~\cite{Fujikawa:1983az} is a pioneering work
on the trace anomaly in lattice gauge theory.}
\begin{align}
   \left\langle
   \mathcal{T}_{\mu\mu}(x)\mathcal{O}
   \right\rangle
   =
   \frac{1}{8}
   \left\langle
   \frac{\partial}{\partial\xi_\alpha}
   \left[\mathcal{Z}\Bar{\Delta}_\xi\gamma_\mu\mathcal{S}_\mu(x)\right]_\alpha
   \mathcal{O}
   \right\rangle
   +8c(N_c^2-1)a^{-4}
   \left\langle
   \mathcal{O}
   \right\rangle,
\label{eq:(3.19)}
\end{align}
where we have used Eq.~\eqref{eq:(3.7)} assuming $x\notin\supp(\mathcal{O})$.

Now we discuss how our energy-momentum tensor~\eqref{eq:(3.4)} can be used in
actual non-perturbative Monte Carlo simulations. As Eqs.~\eqref{eq:(3.3)}
and~\eqref{eq:(3.4)} show, the definition of the energy-momentum tensor
contains four combinations of renormalization constants,
$\mathcal{Z}^2\mathcal{Z}_S$, $\mathcal{Z}_{\text{EOM}}$,
$\mathcal{Z}_T/\mathcal{Z}_S$ and~$c$; we have to know these numbers to
construct the energy-momentum tensor. Among these, the
ratio~$\mathcal{Z}_T/\mathcal{Z}_S$ has been measured~\cite{Donini:1997hh,%
Farchioni:2000mp,Farchioni:2001yn,Farchioni:2004fy,Demmouche:2010sf} by
determining the coefficients in the relation ($\mathcal{O}$ is a fermionic
spinorial operator)
\begin{equation}
   \partial_\mu^S\left\langle S_\mu(x)\mathcal{O}\right\rangle
   +\frac{\mathcal{Z}_T}{\mathcal{Z}_S}
   \partial_\mu^S\left\langle T_\mu(x)\mathcal{O}\right\rangle
   =\frac{1}{\mathcal{Z}_S}\left(M-\frac{1}{a}\mathcal{Z}_\chi\right)
   \left\langle\chi(x)\mathcal{O}\right\rangle.
\label{eq:(3.20)}
\end{equation}
This is Eq.~\eqref{eq:(2.21)} up to~$O(a)$ corrections. Thus we already know
that the number~$\mathcal{Z}_T/\mathcal{Z}_S$ can be determined numerically (up
to~$O(a)$ corrections). Next, to determine~$\mathcal{Z}_{\text{EOM}}$
in~Eq.~\eqref{eq:(3.3)}, we may use the conservation law
of~$\varTheta_{\mu\nu}(x)$, Eq.~\eqref{eq:(3.15)}, itself.

The determinations of~$\mathcal{Z}^2\mathcal{Z}_S$ and~$c$ are somewhat
correlated. We first note that, when the energy-momentum
tensor~$\mathcal{T}_{\mu\nu}(x)$~\eqref{eq:(3.4)} is inserted in a
physical amplitude, the last term of~Eq.~\eqref{eq:(3.4)} just gives rise to an
additive constant~\eqref{eq:(3.7)} that is identical for \emph{any\/} physical
amplitudes; the constant~$c$ thus corresponds to a choice of the origin of the
energy. If we consider the difference in expectation values
of~$\mathcal{T}_{\mu\nu}(x)$ in two physical states, therefore, the contribution
of the last term of~Eq.~\eqref{eq:(3.4)} cancels out and it is the same as the
difference in expectation values
of~$(1/2)[\varTheta_{\mu\nu}(x)+\varTheta_{\nu\mu}(x)]$. Noting this fact, one
may determine $\mathcal{Z}^2\mathcal{Z}_S$, the absolute normalization
of~$\varTheta_{\mu\nu}(x)$ in~Eq.~\eqref{eq:(3.3)}, by setting the difference of
expectation values of the ``energy operator''
$-a^3\sum_{\Vec{x}}\varTheta_{00}(x)$ in two different physical states
to a certain prescribed value.

Although the value of~$c$ influences neither on the conservation law nor on the
difference in expectation values of the energy-momentum tensor, there exists a
physically natural choice of~$c$ in the present supersymmetric system. That
is, we may require that the expectation value of the energy density to vanish
\begin{equation}
   \left\langle\mathcal{T}_{00}(x)\right\rangle=0,
\label{eq:(3.21)}
\end{equation}
when \emph{periodic boundary conditions} are imposed on all the fields. This
requirement fixes
\begin{equation}
   c=-\frac{a^4}{2(N_c^2-1)}
   \left\langle\varTheta_{00}(x)\right\rangle_{\text{periodic boundary conditions}},
\label{eq:(3.22)}
\end{equation}
because of~Eq.~\eqref{eq:(3.7)}. Eq.~\eqref{eq:(3.21)} states that the
derivative of the supersymmetric partition function (i.e., the Witten
index~\cite{Witten:1982df}) with respect to the temporal size of the system
vanishes; this holds only if a choice of the origin of the energy is consistent
with the SUSY algebra. In other words, Eq.~\eqref{eq:(3.21)} is required for
the spatial integral~$-a^3\sum_{\Vec{x}}\mathcal{T}_{00}(x)$ to be the energy
operator appearing in the right-hand side of the SUSY algebra. It is thus a
natural requirement from the perspective of SUSY.\footnote{We emphasize that
Eq.~\eqref{eq:(3.21)} should hold even in a theory in which SUSY is
spontaneously broken (recall that the Witten index is independent of the
temporal size of the system even in such a case). This same idea was adopted
to find an expression of the energy density in a lattice formulation of the
two-dimensional (2D) $\mathcal{N}=(2,2)$
SYM~\cite{Kanamori:2007ye,Kanamori:2007yx,Kadoh:2009rw}. For this 2D system,
fortunately, there exist lattice formulations that possess one exact
fermionic symmetry~\cite{Cohen:2003xe,Sugino:2003yb,Sugino:2004qd,%
D'Adda:2005zk} and, employing this symmetry, one can define the energy density
operator that shows the zero-point energy consistent with SUSY. This energy
density operator on the lattice has been used~\cite{Kanamori:2009dk} to measure
the vacuum energy density associated with a possible dynamical SUSY breaking in
the above 2D system.} This completes our non-perturbative construction of a
symmetric energy-momentum tensor. It is conserved in the quantum continuum
limit and, by choosing~$c$ as~Eq.~\eqref{eq:(3.22)}, the origin of the energy
is consistent with SUSY.

A remaining important issue is whether our symmetric energy-momentum
tensor~\eqref{eq:(3.4)} generates, through a renormalized WT relation, a
correctly-normalized translation on renormalized fields. Although we do not
go into this question in the present paper, assuming the existence of a
renormalized translational invariant theory and considering the uniqueness of
the conserved symmetric energy-momentum tensor in the continuum
limit~\cite{Caracciolo:1989pt}, we believe that the answer is affirmative:
Of course, further consideration is required on this point.

\section{Conclusion}

We have presented a non-perturbative construction of a symmetric
energy-momentum tensor in the lattice formulation of 4D $\mathcal{N}=1$ SYM.
Inspired by the relation~\eqref{eq:(1.9)} in the FZ supermultiplet, we defined
the energy-momentum tensor by a renormalized modified super transformation of a
renormalized SUSY current that generates a finite super transformation on
renormalized operators (under the tuning of the gluino mass). Then, it can be
shown that the lattice energy-momentum tensor is conserved in the quantum
continuum limit. The resulting energy-momentum tensor may be used to measure
physical quantities related with the energy-momentum tensor (such as the
viscosity) by numerical simulations.

A question naturally arises is then what will be resulted from the
consideration on the another relation in the FZ multiplet,
Eq.~\eqref{eq:(1.8)}, a relation between the $U(1)_A$ current and the SUSY
current. It is clear that the classical relation~\eqref{eq:(1.8)} as it stands
cannot hold in quantum theory, because the $U(1)_A$ current is not conserved
owing to the axial anomaly, while we are believing that SUSY does not suffer
from any anomaly. Ironically, it turns out that the situation for this
seemingly simple relation~\eqref{eq:(1.8)} is much more complicated than for
(seemingly complicated) Eq.~\eqref{eq:(1.9)}; this is precisely because of the
axial anomaly. We hope to come back this problem in the near future. If the
full structure of the FZ multiplet (with corrections by the anomaly) can be
realized in a well-regularized framework such as the lattice, it will be a
quite useful starting point to understand the so-called anomaly
puzzle~\cite{Clark:1978jx,Grisaru:1978vx,Grisaru:1985yk,Grisaru:1985ik,%
Shifman:1986zi,Ensign:1987wy,ArkaniHamed:1997mj,Huang:2010tn,%
Komargodski:2010rb,Yonekura:2010mc,Yonekura:2012uk}.

\section*{Note added in proof}

Eq.~\eqref{eq:(3.13)} implies the following very suggestive relation: Noting
that $\sum_{y\in\mathcal{D}_x}\partial_\nu^S\mathcal{S}_\nu(y)$ in the first term
of the right-hand side of~Eq.~\eqref{eq:(3.13)} does not have the support
at~$y=x$, using the SUSY WT relation~\eqref{eq:(2.25)} again, we have
\begin{align}
   &a^4\sum_{x\in\mathcal{D}_{\mathcal{O}}}
   \left\langle
   \partial_\mu^S\varTheta_{\mu\nu}(x)
   \mathcal{O}\right\rangle
\notag\\
   &=-\frac{1}{8}
   (C^{-1}\gamma_\nu)_{\alpha\beta}\,
   a^4\sum_{x\in\mathcal{D}_{\mathcal{O}}}
   a^4\sum_{y\in\mathcal{D}_x}
\notag\\
   &\qquad{}\times
   \left\langle
   \left\{
   \mathcal{Z}\left[
   -\frac{1}{a^4}\frac{\partial}{\partial\Bar{\xi}(x)}\Delta_\xi
   +a\mathcal{E}(x)\right]_\alpha,
   \mathcal{Z}\left[
   -\frac{1}{a^4}\frac{\partial}{\partial\Bar{\xi}(y)}
   \Delta_\xi+a\mathcal{E}(y)\right]_\beta\right\}
   \mathcal{O}\right\rangle,
\end{align}
where $\mathcal{D}_{\mathcal{O}}$ is a finite region that contains the
operator~$\mathcal{O}$ entirely. On the other hand, for the anti-symmetric part
of~$\varTheta_{\mu\nu}(x)$, $\mathcal{A}_{\mu\nu}(x)$ in~Eq.~\eqref{eq:(3.16)},
\begin{equation}
   a^4\sum_{x\in\mathcal{D}_{\mathcal{O}}}
   \left\langle
   \partial_\mu^S\mathcal{A}_{\mu\nu}(x)
   \mathcal{O}\right\rangle\xrightarrow{a\to0}0,
\end{equation}
because Eq.~\eqref{eq:(3.17)} shows that
$\langle\mathcal{A}_{\mu\nu}(x)\mathcal{O}\rangle$ is proportional to the
delta function~$\delta^4(x-z)$ in the continuum limit, where $z$ is any point
in the support of~$\mathcal{O}$. Similarly, for the last term
of~Eq.~\eqref{eq:(3.4)}, we have
\begin{equation}
   a^4\sum_{x\in\mathcal{D}_{\mathcal{O}}}
   \left\langle
   \partial_\mu^S\delta_{\mu\nu}
   \tr\left[\Bar{\psi}(x)(D+M)\psi(x)\right]
   \mathcal{O}\right\rangle=0.
\end{equation}
Thus, combining the above three relations,
\begin{align}
   &a^4\sum_{x\in\mathcal{D}_{\mathcal{O}}}
   \left\langle
   \partial_\mu^S\mathcal{T}_{\mu\nu}(x)
   \mathcal{O}\right\rangle
\notag\\
   &\xrightarrow{a\to0}-\frac{1}{8}
   (C^{-1}\gamma_\nu)_{\alpha\beta}\,
   a^4\sum_{x\in\mathcal{D}_{\mathcal{O}}}
   a^4\sum_{y\in\mathcal{D}_x}
\notag\\
   &\qquad{}\times
   \left\langle
   \left\{
   \mathcal{Z}\left[
   -\frac{1}{a^4}\frac{\partial}{\partial\Bar{\xi}(x)}\Delta_\xi
   +a\mathcal{E}(x)\right]_\alpha,
   \mathcal{Z}\left[
   -\frac{1}{a^4}\frac{\partial}{\partial\Bar{\xi}(y)}
   \Delta_\xi+a\mathcal{E}(y)\right]_\beta\right\}
   \mathcal{O}\right\rangle.
\end{align}
This may be regarded as the SUSY algebra in the present lattice framework if
we assume that the lattice energy-momentum tensor~$\mathcal{T}_{\mu\nu}$
generates the translation (on gauge invariant operators~$\mathcal{O}$) as
\begin{align}
   a^4\sum_{x\in\mathcal{D}_{\mathcal{O}}}
   \left\langle
   \partial_\mu^S\mathcal{T}_{\mu\nu}(x)
   \mathcal{O}\right\rangle
   \xrightarrow{a\to0}
   -\left\langle\partial_\nu\mathcal{O}\right\rangle.
\end{align}

\section*{Remark in proof}
A quite similar but somewhat different definition of a lattice energy-momentum
tensor in 4D $\mathcal{N}=1$ SYM, which is superior in several aspects compared
with the one in the present note, has been given in~Ref.~\cite{Suzuki:2012wx}.

\section*{Acknowledgements}

I would like to thank Michael G.~Endres, Kazuo Fujikawa, Yoshio Kikukawa,
Martin L\"uscher, Yusuke Taniguchi and Naoto Yokoi for enlightening
discussions. Some part of this work was carried out during the YITP workshop on
``Field Theory and String Theory'' (YITP-W-12-05) and I would like to thank the
Yukawa Institute for Theoretical Physics at Kyoto University for the
hospitality. This work is supported in part by a Grant-in-Aid for Scientific
Research, nos.~22340069 and~23540330.

\bibliographystyle{elsarticle-num}


\section{Introduction}

When confined to two-dimensional (2D) space, condensed matter behaves
differently than in three dimensions (3D). A striking
example is the physics of the Kosterlitz-Thouless transition in superfluid
films, not to mention the phenomenon of high-temperature superconductivity
in the cuprates, which is intimately
related to the motion of electrons within weakly-coupled copper-oxide layers.
Another example is provided by 2D crystals, where thermal fluctuations are
so strong as to rule out long-range translational (but not orientational)
order for non-zero temperatures, leaving it open for the possibility of
unconventional melting scenarios. The celebrated KTHNY theory of 2D melting
thereby predicts a continuous two-stage melting from
a crystalline to a hexatic
phase and, subsequently, from the hexatic to an isotropic liquid~\cite{KTHNY}.
The intermediate hexatic phase has short-range translational order
(i.e., it is fluid-like) but quasi-long-range orientational order,
characterized by a power-law decay of bond-angle pair correlations.
Other possibilities are a standard first-order 2D-melting transition
(as in Chui's theory~\cite{Chui}) or a stronger, i.e., first-order
hexatic-to-isotropic transition, as shown by Bernard and Krauth
to be the case for hard discs~\cite{Bernard}. In some cases, melting is
one-stage and first-order but the hexatic phase can be accessed via a
metastable route~\cite{Chen,Mejia}.

The KTHNY theory has been confirmed many times in simulation and experiment,
especially for particles with long-ranged interactions, though in a few cases
with some controversy about the order of the two transition
steps~\cite{Murray,Marcus,Zahn,Keim,Lin1,Han,Peng,Muto,Lin2,Lee,Qi,Clark}.
From a computational point of view, the assessment of the
order of 2D melting in practical cases can be hard, due to important
finite-size effects and slow relaxation to equilibrium. An extreme case
is hard discs where the exact nature of the melting transition could be
established only recently, by employing unprecedentedly large system sizes
of the order of a million particles in a box of fixed
volume~\cite{Bernard}. Generally speaking, working in the isothermal-isobaric
ensemble (rather than at constant volume) can ease the interpretation
of the simulation data, since one cannot incur in any of the
finite-size artifacts that harm constant-volume simulations.
Recently, we have provided unambiguous evidence of the occurrence of
continuous melting via a hexatic phase in the 2D Gaussian-core model
(GCM)~\cite{Prestipino1}. This system, taken as a prototype of the
phase behavior of ``empty'' polymers~\cite{Likos1}, has long been known
to exhibit reentrant melting (i.e., solid melting under isothermal
compression) as well as water-like anomalies in three
dimensions~\cite{Stillinger,Lang,Prestipino2,Mausbach}. In two dimensions,
the GCM melting becomes two-staged, with an extremely narrow hexatic
region whose properties comply with the predictions of the KTHNY theory.
In particular, besides a reentrant isotropic fluid, also a reentrant
hexatic phase exists. Recently, water-like anomalies have been reported
even in the one-dimensional GCM~\cite{Speranza}, where however no properly
defined thermodynamic transition occurs.

A class of systems which exhibit reentrant melting, solid polymorphism,
and other water-like anomalies is formed by particles interacting via
spherically-symmetric core-softened (CS) potentials~\cite{Buldyrev}.
For these systems, the strength of the repulsive component of the
interaction undergoes some sort of weakening over a range of interparticle
distances. The effects of particle-core softening on thermodynamic behavior
were originally investigated by Hemmer and Stell~\cite{Hemmer}, who were
interested in the possibility of multiple critical points and
isostructural solid-solid transitions in simple-fluid systems.
A few years later, Young and Alder~\cite{Young} showed that the phase diagram
of the hard-core plus square-shoulder repulsion exhibits an anomalous melting
line of the kind observed in Cs and Ce. Later, Debenedetti~\cite{Debenedetti}
showed that systems of particles interacting via continuous CS potentials
are capable of contracting when heated isobarically (a behavior which goes
under the name of ``density anomaly'').
In the last decade or so, there has been a renewal of interest in the phase
behavior of CS potentials, which has led to the discovery of many unusual
properties such as multiple reentrant-melting lines, polymorphism both in
the fluid and solid phases, stable cluster solids (for bounded interactions
only), and a plethora of thermodynamic and dynamic water-like
anomalies~\cite{Sadr,Jagla,Likos2,Kumar,Deoliveira,Gibson,Lomba,Fomin,Xu,Pizio,Pamies,Pauschenwein,Malescio1}.
A common feature of CS potentials is their ability to generate two distinct
length scales in the system, one related to the inner core and the other
associated with the milder component of the repulsion~\cite{Buldyrev}.
Due to this, CS fluids are characterized by two competing, expanded and
compact, local particle arrangements. Such an interplay has a disruptive
influence on crystal order, the fluid phase being thus recovered upon
compression.
This two-scale or two-fluid property mimics the behavior of the more complex
network-forming fluids, where the looser and denser local structures arise
from the incessant building and breaking of the dynamic network generated
by directional bonds. As a result, a paradigm has been established whereby
the existence of two competing local structures is essential for the
occurrence of anomalous behaviors in simple fluids.

Unusual phase behavior has also been found in 2D systems.
Scala and coworkers~\cite{Scala} carried out molecular-dynamics simulations
of the 2D square-shoulder plus
square-well (SSSW) potential finding water-like anomalies and identifying
two different solid phases: a triangular crystal (low-density phase) and a
square crystal (high-density phase).
However, the observation of a hysteresis loop in their simulations
suggested that the liquid-solid transition is actually first-order,
thus washing out the possibility of a continuous melting to a hexatic phase.
Wilding and Magee also studied the 2D SSSW potential and showed that the
thermodynamic
anomalies of the model, rather than stemming from a metastable liquid-liquid
critical point as previously surmised, were induced by the quasicontinuous
nature of the 2D freezing transition~\cite{Wilding}.
Almudallal and coworkers~\cite{Almudallal} contested this result, showing by
various free-energy techniques that all transitions in the 2D SSSW model are
in fact first-order.
Malescio and Pellicane studied a 2D system of particles interacting
through a potential consisting of an impenetrable hard core plus a square
shoulder. They found a variety of stripe phases whose formation was imputed
to the existence of two characteristic length scales~\cite{Malescio2}.
All such investigations suggest that the existence of two length scales
would be a prerequisite for observing water-like anomalies also in 2D.

In a series of recent papers the two-fluid
picture has been overthrown and a novel minimal scenario for the occurrence of
anomalies has been established~\cite{Saija,Prestipino3,Malescio3,Malescio4}.
More specifically,
it has been shown that a weakly-softened isotropic pair repulsion, with
a single characteristic length which becomes more loosely defined in a range
of pressures, is able to give rise to an unusual phase behavior.
Our purpose here is to verify whether anything similar occurs in 2D, where
in addition an interesting interplay with hexatic order may take place.
Far from being purely academic, the present study can be relevant for many
soft materials. For instance, one monolayer of $N$-isopropylacrylamide
(NIPA) microgel spheres confined between two glass cover slips is an ideal
system to study 2D melting because the effective interparticle potential
is short-ranged and repulsive, with a temperature-tunable volume
fraction~\cite{Dullens,Lin2,Han,Peng,Alsayed}.
Systems like this would be natural candidates for detecting, by video optical
microscopy or light scattering, the kind of phenomena that we are going to
illustrate below.

The plan of the paper is the following. After introducing our 2D model and
method in Section II, we sketch the model phase diagram in Section III,
highlighting the unconventional structure of the fluid phase which, similarly
to other instances of CS repulsion, shows a number of water-like anomalies
without any interplay between two characteristic length scales.
The melting mechanism is studied in more detail in Section IV,
where we show the existence of a hexatic phase.
Section V is finally devoted to concluding remarks.

\section{Model and method}
\setcounter{equation}{0}
\renewcommand{\theequation}{2.\arabic{equation}}

We consider a purely repulsive pair potential in two dimensions, modelled
through an exponential form which was first introduced, about four decades ago,
by Yoshida and Kamakura (YK)~\cite{Yoshida}:
\begin{equation}
u(r)=\epsilon\exp\left\{a\left(1-\frac{r}{\sigma}\right)-
6\left(1-\frac{r}{\sigma}\right)^2\ln\left(\frac{r}{\sigma}\right)\right\}\,,
\label{2-1}
\end{equation}
where $\epsilon$ and $\sigma$ set the energy and length scales, respectively,
and $a\ge 0$.
The YK potential behaves as $r^{-6}$ for small $r$, and falls off very rapidly
for large $r$. The smaller $a$ the higher the degree of softness of
$u(r)$, i.e., the flatter the repulsive ``shoulder'' around $r=\sigma$.
Technically speaking~\cite{Debenedetti}, $u(r)$ can be regarded as soft only
for $a\lesssim 2.3$, since only in this range there is an
interval of distances where the local virial function $-ru'(r)$ decreases
when the interparticle separation decreases.
However, it is over the wider range $a\lesssim 5.5$ that the phase diagram
is expected to show a reentrant-fluid region~\cite{Prestipino3,Prestipino4}.
We shall here focus on the case $a=3.3$ (Fig.\,1), which was already shown to
possess a rich anomalous phase behavior in three dimensions which cannot
simply be explained by the existence of two distinct nearest-neighbor (NN)
distances in the dense fluid~\cite{Prestipino3}.

We now briefly describe the methodology of the present investigation.
After identifying the relevant crystal phases by means of total-energy
calculations at zero temperature
(of the type illustrated in, e.g., Ref.\,\cite{Prestipino5}),
we explore the phase diagram by Monte Carlo (MC) simulation in the $NPT$
ensemble. As is common practice, we adopt periodic boundary conditions and
employ cell linked lists in order to speed up the simulation.
While systems of about 1000 particles are suitable for investigating
bulk properties and for determining the approximate location of
phase boundaries, we consider larger systems of about 6000 particles
for the search and characterization of the hexatic phase at selected
pressures. For the same $P$ values, we check the order of the melting
transition independently through thermodynamic integration combined
with exact free-energy calculations.
For the fluid phase, a dilute gas is used as a reference state,
whose chemical potential is calculated by the Widom method, whereas a
low-temperature crystal is chosen as the starting point of the simulation in
the solid region of the phase diagram. At this initial state, the excess
Helmholtz free energy of the system is computed by the Einstein-crystal
method~\cite{Frenkel}.

\section{Fluid structure and water-like anomalies}

For $a=3.3$, the repulsive shoulder at $\sigma$, which is still visible
in the plot of $u(r)$ for $a\simeq 2$, has by then completely faded out,
surviving only in the form of a modest bump in the otherwise
monotonously-decreasing profile of the virial function (see Fig.\,1).
Yet, this almost structureless potential shows three distinct stable crystal
arrangements at $T=0$ in two dimensions, as it follows from the
calculation of the chemical potential $\mu$ as a function of the pressure $P$
for all five Bravais lattices and for the honeycomb lattice.
Upon increasing $P$, the sequence of phases is triangular-square-triangular,
the square crystal being thermodynamically stable in the pressure range
$2.27$-$4.36$, in units of $\epsilon/\sigma^3$ (from now on, pressure and
temperature will be given in reduced units).

Assuming that no other crystal phases come into play at nonzero temperatures,
we first sketched the overall phase diagram by the heat-until-it-melts
(HUIM) method. In practice, for selected $P$ values we run a chain of MC
simulations of the solid stable at the given $P$, for increasing $T$ values
at regular intervals of 0.0005, until we observe a clear
jump in both the particle-number density $\rho=N/V$ and total energy per
particle $E/N$, which signal the melting of the solid. The character of this
transition will be addressed in the next section. The HUIM method simply
overlooks the possibility of solid overheating, but this is usually a
reasonable
approximation if one is not interested in very precise estimates of the
melting temperature $T_m$ ({\it a posteriori}, the typical
error implied by the HUIM method in determining the melting point of the
present system was about 10\%, and smaller for the square crystal).
The outcome of this analysis
is reported in Fig.\,2 where the most revealing feature is the reentrance
of the fluid phase as the system, already settled in the low-density
triangular solid, is further compressed at not too low temperatures.
Fluid reentrance occurs even twice in a smaller range of temperatures above
0.020.

As we already said in the Introduction, until very recently the conventional
wisdom on the origin of the reentrant-melting phenomenon in CS systems
with an unbounded interparticle repulsion rested on a competition,
with destabilizing effects on crystal ordering, between two ways of
arranging particles close to each other, which gives rise to either an
expanded or a more compact structure in the dense fluid. This may only occur
provided the pressure is strong enough as to bring neighboring particles at
distances close to $\sigma$. Clearly, this explanation cannot work for
a case like the present one, where it is hard to maintain that there are
two different length scales in the potential. Rather, we may view the
question from a solid-state perspective and argue that a succession of
reentrant-melting lines in the phase diagram is simply the outcome of
the existence of multiple solid phases at low temperature. On increasing
pressure within the range of stability of any of such solids, the crystal
strength first grows up to a maximum (and in parallel also $T_m$),
but then it progressively reduces on approaching the boundary
of the next stable solid at higher pressure. After all, this is exactly
the rationale behind the melting criterion discussed in
Ref.\,\cite{Prestipino4},
which in fact is very effective for soft repulsive interactions.

We found confirmation that the present model shows only one repulsive
length scale by computing the fluid radial distribution function (RDF)
along the isotherm at $T=0.1$, i.e., just above the maximum $T_m$
of the low-density triangular solid (see Fig.\,3). Looking
at the position of the main RDF peak, we see a systematic shift to lower
and lower distances upon compression, with a slight widening of the peak
around $P\simeq 2$, i.e., next to the first reentrant-melting line.
We interpret this as the evidence of a unique NN characteristic distance
in the system, which becomes more loosely defined across the region of
reentrant melting. This should be contrasted with what occurs for a more
conventional CS repulsion, where the first RDF maximum is twin-peaked,
with two definite NN distances that take turns at providing the absolute
maximum for the RDF.

Aside from the existence of one repulsive length scale rather than two,
many water-like anomalies are also found in the present
system, starting with a line of $\rho$ maxima in the
fluid close to the first reentrant-melting line (see Fig.\,2).
This is a typical occurrence in systems with CS potentials,
where however the region of volumetric anomaly is usually more extended.
We checked that the line of the density anomaly ceases to exist just before
plunging into the solid region. Another type of anomaly is the
so-called structural anomaly, that is a non-monotonous pressure behavior
of the amount of ``translational order'' in the fluid, as measured through
the value of minus the pair entropy per particle ($-s_2$)~\cite{Prestipino6}
(see Fig.\,2 and the left top panel of Fig.\,4).
Rather than monotonously increasing with pressure at constant temperature,
the degree of spatial order reduces in the reentrant-fluid region, as a
result of a looser definition of the NN distance.
A non-sharp average separation between neighboring particles
acts as a perturbing factor for the local order, bringing about a slight
decrease of $-s_2$ with pressure.
In the same range of pressures where $-s_2$ decreases,
the self-diffusion coefficient $D$, which we
measure by a series of $TVN$ molecular-dynamics runs, gets enhanced with
pressure (Fig.\,2 and right top panel of Fig.\,4). Both lines of structural
and diffusional anomalies appear to sprout out of the point of maximum $T_m$
for the low-density triangular solid and roughly terminate at the point of
maximum $T_m$ for the square solid. This is again similar
to other CS systems, also for what concerns the crossing of the anomaly lines
(see e.g. Fig.\,11 of \cite{Jabes}). Finally, we checked for just one pressure
value ($P=2$) the existence of a temperature minimum in both the isobaric
specific heat $C_P$ and isothermal compressibility $K_T$ (lower panels
of Fig.\,4), similarly to what found for water at ambient pressure.
Both minima fall in the high-temperature fluid region,
much far from the maximum-density point for the same pressure.

\section{Hexatic behavior}

The existence in our model of two solids with different crystal symmetry,
i.e., triangular and square, gives the opportunity to investigate the
intriguing
possibility of two distinct intermediate phases (hexatic and tetratic)
between the solid and the normal fluid, with the further bonus of a
yet-to-be-observed
hexatic-tetratic transition near the confluence point between the two solids
and the fluid.

In order to clarify the melting scenario for a given $P$,
we run our computer code along two simulation paths, one starting
from the solid phase at $T=0.005$ and the other from a high-temperature
fluid state. We advance in steps of $\Delta T=0.005$, equilibrating the
system for long before generating an equilibrium trajectory of $M$ sweeps
(one sweep corresponding to $N$ trial MC moves), with $M$ ranging between
$5\times 10^5$ and $3\times 10^6$, depending on how far we are from melting.
Once a guess of the transition point is made, we restart the simulation
slightly ahead of it with a smaller $\Delta T$ and/or a larger $M$
in order to better discriminate between first-order and continuous melting.
Besides $\rho$ and $E/N$, we measure two order parameters (OP),
$\psi_T$ and $\psi_O$, which are sensitive to the overall translational
and orientational triangular/square order, respectively. The precise
definition of both quantities has been given in Ref.\,\cite{Prestipino1},
with obvious modifications for the square-lattice case. Moreover, we keep
track of the OP susceptibilities $\chi_T$ and $\chi_O$, defined as $N$
times the variance of the respective OP estimator. Finally, we calculate
two orientational correlation functions (OCF)~\cite{Prestipino1}, $h_6(r)$
and $h_4(r)$, which inform on the typical size of a space region in which
NN-bond angles are strongly correlated. The KTHNY theory predicts an
algebraic $r^{-\eta(T)}$ large-distance decay of the OCF in the hexatic
phase, at variance with what occurs in an isotropic fluid where the decay
is much faster, i.e., exponential. According to the same theory, $\eta$
equals $1/4$ at the transition point between hexatic and isotropic fluid.

In Figs.\,5-8, we report $\rho$ and $E/N$ for two different system sizes
and various simulation protocols, as a function of $T$ for
$P=0.5,2,3$, and 5. We clearly see that, while melting is
continuous for $P=0.5,2$, and 5, it is
certainly first-order for $P=3$ (and 4, data not shown), as evidenced by the
hysteresis loops. Based on our experience with the GCM~\cite{Prestipino1},
we can conjecture that there is a narrow region of hexatic phase for
$P=0.5,2$, and 5
(to be confirmed later by the analysis of OPs and OCF), whereas no
definite conclusion can be reached at this point for $P=3$, where a
tetratic phase can in principle exist even in presence of a first-order
transition. A special remark is due for $P=5$, where,
in contrast with what occurs for smaller pressures, the crystal energy
decreases upon heating, which means that the increase in kinetic energy
is more than compensated for by the loss in potential energy, whose high
rate of decrease is due to a NN distance lying in the harsh part of the
potential core.

For three pressures ($0.5,2$, and 3),
we plot the OPs and related
susceptibilities in Figs.\,9-11. For the first two pressures, we see that
$\psi_T$ vanishes at a slightly smaller temperature than $\psi_O\equiv\psi_6$,
which implies that the hexatic phase is confined to a narrow $T$ interval
not wider than 0.0005 for $P=0.5$ (0.0015 for $P=2$), as also witnessed
by the maxima of the two susceptibilities occurring at slightly different
$T$ values. These temperature intervals compare well with the $T$ range
of the bridging region between the solid and fluid branches in Figs.\,5
and 6. We thus confirm the same phenomenology of the GCM, namely that the
width of the hexatic region increases with pressure.
For $P=5$, the findings are similar to $P=2$, with roughly
the same width of the hexatic region and comparable levels of orientational
and translational order in the hot solid (data not shown).
Going to $P=3$, the picture is quite different since the two OPs now apparently
vanish at the same temperature; at that point, the orientational
susceptibility shows a spike (rather than a critical peak)
which is usually associated with a first-order transition.
All evidence suggests that standard first-order melting is a plausible
explanation for $P=3$, and the same conclusion can be made for
$P=4$ (data not shown).

A more direct evidence of the existence of a bond-angle ordered fluid
(or a clue to its absence) comes from the large-distance behavior of
the OCF. We plot this function in Figs.\,12-14
for various temperatures across the relevant
region, for $P=0.5,2$, and 3, respectively. It appears that, for $P=0.5$
and 2, $h_6(r)$ decays algebraically in a $T$ region of limited extent,
roughly corresponding to the bridging region in Figs.\,5 and 6.
Moreover, the decay exponent in the hexatic region is smaller than
$1/4$, becoming larger only on transforming to the isotropic fluid,
and the same is found for $P=5$ (data not shown).
On the contrary, for $P=3$ and 4, $h_4(r)$ switches directly from no decay
at all to an exponential damping, showing that there is no tetratic phase
in our model.

The obvious question now arises as to why, at variance with the hexatic one,
the tetratic phase is not stable in our model. This question is difficult
to answer since
it involves the consideration of the delicate equilibrium between the energy
and entropy of the two competing phases, in this case the tetratic phase and
the isotropic fluid phase. A possible hint could be the level of orientational
order that is found in the two types of solid slightly before melting.
If we look at Figs.\,9-11, we see that the value of $\psi_O\equiv\psi_4$
for $P=3$ is less than a
half of $\psi_6$ (note that the same holds for $P=4$), and
this would explain why, in the square-crystal case, long-range orientational
order does not survive (even in the weakened form typical of a tetratic
phase) the loss of quasi-long-range positional order determined by melting.

Finally, we checked by an indipendent route the order of the melting
transition for $P=0.5,2,3$, and 4. As anticipated, for each pressure we
carried out MC simulations along two different paths, one beginning
from a cold solid and the other from a low-density fluid.
Using thermodynamic integration in combination with exact free-energy
calculations at the initial points of the paths, we were able to obtain
the system chemical potential $\mu$ along the solid and fluid branches
as a function of $T$. For $P=0.5$ and 2, we did not find any crossing
of the two $\mu(T)$ curves, thus confirming a continous melting transition.
Notwithstanding the care we put in keeping under control any source of
statistical error, we nevertheless found a small discrepancy (about
$0.0005\epsilon$, practically constant in a 0.01 wide $T$ range
across the melting transition) between the $\mu$ levels of the solid
and the fluid, which is presumably due to the not-so-small $\Delta T$
employed along the paths far from melting. On the contrary, for $P=3$
and 4, we observed a clear crossing between the two $\mu(T)$ curves,
respectively at $T=0.0269$ and $T=0.0176$, which is consistent with
the location of the density and energy jumps (for $P=3$, see Fig.\,7),
and suggestive of a first-order transition.

Summing up, we collected multiple evidence of a continuous
melting via a hexatic phase for the triangular crystals,
while the melting transition is certainly discontinuous and ``standard''
for the square crystal. Nothwithstanding the ``small''
sizes of the investigated samples, we think that our conclusions are robust
since they result from many independent indicators of the phase-transition
order.

\section{Conclusions}

We have analyzed the phase behavior of purely-repulsive 2D particles
with a weakly-softened core. The nature of this repulsion is such as to
determine a characteristic nearest-neighbor distance in the fluid phase
whose statistical precision, expressed by the width of the main peak of
the radial distribution function, shows a non-monotonous behavior with
pressure at not too high temperatures. Notwithstanding the fact that the
two-fluid paradigm of core-softened (CS) potentials does not apply here,
we anyway observe the same phenomenology as in conventional CS systems,
with solid polymorphism, multiple reentrant-melting lines, and many other
water-like anomalies.
While this is similar to the 3D case~\cite{Prestipino3}, the melting
transition of the present system is different, since it is continuous
and occurs via an (even reentrant) hexatic phase for the triangular
solids while being standard first-order for the
intermediate square solid (i.e., no parallel tetratic
phase exists). The hexatic behavior appears to be
consistent with the KTHNY theory,
as witnessed by the value of the decay exponent of the orientational
correlation function at the hexatic-to-isotropic fluid transition temperature.

Our findings could be relevant for many real soft-matter systems.
Already nowadays, monodisperse colloidal suspensions can be engineered
in such a way as to exhibit a temperature-modulated repulsion with some
amount of softness. Probably, the day is not far when by appropriate
functionalization of the particle surface a colloidal system will be shown
to exhibit the kind of phase-transition phenomena that we have illustrated
here for a rather specific model potential.