\section{Introduction}

It has long been realized that deeply inelastic scattering (DIS) off polarized protons probes the physics of the chiral anomaly in QCD~\cite{Kodaira:1979pa} though its precise role has been the subject of some debate~\cite{Altarelli:1988nr,Carlitz:1988ab,Jaffe:1989jz}. The purpose of this work is revisit and cast new light on the chiral anomaly with a view to better understand the interplay between parton dynamics and the topology of the QCD vacuum in the helicity structure of the proton at high energies.

In this paper, we will focus on the triangle graph~\cite{Adler:1969gk,Adler:1969er,Bell:1969ts,Kogut:1974kt} whereby the anomaly manifests itself in the coupling of the isosinglet axial vector current to the topological charge density in the polarized proton. A careful treatment of the triangle graph is essential to a first principles understanding of polarized DIS. For instance, as we shall discuss, the off-forward matrix element for the polarized $g_1$ structure function contains an infrared pole that appears to diverge in the forward limit~\cite{Jaffe:1989jz}. It is well known that the triangle graph is embedded in the usual box diagram for polarized DIS in the Bjorken limit of large squared momentum transfer $Q^2$. Our analysis, performed in a worldline formalism  particularly suited to discussions of the anomaly~\cite{Polyakov:1987ez,AlvarezGaume:1983ig,Strassler:1992zr,DHoker:1995aat,DHoker:1995uyv,Haack:1998uy,McKeon:1998et,Schubert:2001he,Bastianelli:2003bg}, will show however that the usual operator product expansion (OPE) formalism is not necessary for this result though it is sufficient. 

Our novel result is that the triangle graph appears identically in the box diagram for high energy polarized DIS in the $x_{B}\ll 1$ Regge asymptotics~\cite{Jaffe:1996zw} of the Bjorken variable $x_{B}$. While there have been  qualitative discussions~\cite{Bass:1997zz,Bass:2004xa,Wakamatsu:2019ain} of the triangle anomaly at small $x_{B}$, a quantitative discussion has been lacking thus far. Our aim is to redress this lack especially in view  of polarized DIS experiments  at the Electron Ion Collider (EIC) in the near future that will access very small values of $x_{\rm B}$ for the first time~\cite{Accardi:2012qut,Aschenauer:2017jsk,Aschenauer:2020pdk}. 

Albeit our focus here is on the triangle anomaly, follow-up papers (Papers II \& III) in preparation will discuss further the fundamental issues underlying the non-perturbative regularization of the infrared pole of the anomaly. Some of these issues were discussed previously 
 by Shore and Veneziano~\cite{Shore:1990zu,Shore:1991dv}, and by Shore, Narison and Veneziano~\cite{Narison:1994hv,Narison:1998aq}--for a nice review, see \cite{Shore:2007yn}. Specifically, in Paper II, we will motivate in the worldline formalism an effective action for Regge asymptotics that is  consistent with anomalous chiral Ward identities~\cite{Wess:1971yu}. In Paper III, we will  discuss the energy evolution of helicity dependent distributions in this framework. We note that there is a considerable body of work on perturbative resummations of the large logarithms in $x_B$ that drive the energy evolution of helicity dependent distributions in Regge asymptotics~\cite{Kirschner:1983di,Bartels:1995iu,Bartels:1996wc,Kovchegov:2015pbl,Kovchegov:2016weo,Kovchegov:2017jxc,Chirilli:2018kkw,Boussarie:2019icw,Cougoulic:2019aja,Kovchegov:2020hgb,Cougoulic:2020tbc}.


To proceed further, we will recap briefly the discussion of the chiral anomaly in polarized DIS and some of the subtle issues in its interpretation. Polarized inclusive deeply inelastic scattering (DIS) is defined as the process 
\begin{eqnarray}
l(l) + N(P, S)\to l(l') + X\,,
\end{eqnarray}
where the lepton $l$ interacts with the polarized target hadron $N$ via the exchange of a virtual photon $\gamma^\ast$ with momentum $q = l - l'$. Here the target is characterized by its momentum vector $P = (P^+, M^2/2P^+, 0_\perp)$ and the covariant spin vector satisfies $S^2 = -1$.

The hadron tensor in DIS is the matrix element of the product of electromagnetic currents~\cite{Blumlein:2012bf},
\begin{eqnarray}
W^{\mu\nu}(q, P, S) = \frac{1}{2\pi} \int d^4x \,e^{iqx} \langle P,S| j^\mu(x) j^\nu(0)|P,S\rangle\,,
\label{hT}
\end{eqnarray}
where $j^\mu = \sum_f e_f \bar{\Psi}_f\gamma^\mu \Psi_f$ is bilinear in the quark and antiquark field operators and  $e_f$ denotes the electric charge of a quark of flavor $f$.
It can be split into symmetric and antisymmetric parts as
\begin{eqnarray}
W^{\mu\nu}(q, P, S) = \bar{W}^{\mu\nu}(q, P) + i \tilde{W}^{\mu\nu}(q, P, S)\,.
\label{SAT}
\end{eqnarray}
Since our interest in this paper is on spin effects in DIS,  our focus will be on the antisymmetric part of Eq.~(\ref{SAT}), which can be expressed in terms of spin dependent structure functions~\cite{Anselmino:1994gn} as
\begin{eqnarray}
\tilde{W}_{\mu\nu}(q, P, S)= \frac{2M_N}{P\cdot q}\epsilon_{\mu\nu\alpha\beta} q^\alpha\Big\{ S^\beta g_1(x_B, Q^2) + \Big[S^\beta - \frac{(S\cdot q)P^\beta}{P\cdot q}\Big]g_2(x_B, Q^2)\Big\}\,,
\label{WA}
\end{eqnarray}
where the virtuality of the incoming virtual photon $Q^2 = -q^2 > 0$,  the Bjorken variable $x_B = Q^2/(2P\cdot q)$, $M_N$ is the proton mass and the totally antisymmetric Levi-Civita tensor $\epsilon_{\mu\nu\alpha\beta}$  is defined with $\epsilon_{0123} = -1$. It is convenient to consider a longitudinally polarized target with the covariant spin vector $S^\mu(\lambda) \simeq \frac{2{\tilde \lambda}_P}{M_N}P^\mu$, where ${\tilde \lambda}_P = \pm \frac{1}{2}$ is the helicity. In this case, the $g_2$ structure function does not contribute.

In the parton model, at leading twist, this expression simplifies to read~\cite{Leader:2001gr}
\begin{equation}
\label{eq:parton-g1}
g_1(x_B, Q^2) = \frac{1}{2} \sum_f e_f^2 \left( \Delta q_f(x_B,Q^2) +\Delta {\bar q}_f(x_B,Q^2)\right)\,,
\end{equation}
where the polarized parton distribution function (pdf)
\begin{equation}
\Delta q_f(x_B,Q^2) = \frac{1}{4\pi} \int dy^- \,e^{-iy^- x_B \,P^+}\,\langle P,S| {\bar \Psi_f}(0,y^-,0_\perp) \gamma^+ \gamma_5 \Psi_f(0)|P,S\rangle\,.
\end{equation}
Here $P^+$ is the large light cone component of the momentum of the hadron. In light front quantization, $\Delta q_f(x_B,Q^2)$ has the physical interpretation of the difference in the density of left and right 
handed quarks of a given quark flavor. Likewise, $\Delta {\bar q}_f(x_B,Q^2)$ denotes the difference in the density of left and right handed anti-quarks of the given flavor. 

The first moment of Eq.~(\ref{eq:parton-g1}) can be expressed~\cite{Ellis:1973kp}, assuming flavor $SU(3)$, as
\begin{equation}
\int_0^1 dx_B\, g_1(x_B,Q^2) = \frac{1}{18} \left(3 F + D +
2 \,\Sigma(Q^2)\right) \,.
\label{g1fm}
\end{equation}
Here $F$ and $D$ in the combinations $F+D$ and $3 F-D$ are proportional respectively to the isotriplet axial vector current and the octet axial vector current. The former is nothing but $g_A$ the 
nucleon's axial vector coupling and is determined quite precisely from $\beta$-decay experiments. Likewise, the latter is well known from hyperon decay experiments. Their running is very weak and they can be treated 
for all relevant purposes as constants. 

The object of interest in this equation is the net light quark helicity  $\Sigma(Q^2)$ defined as the flavor singlet sum:
\begin{equation}
\Sigma(Q^2) = \sum_f \int_0^1 dx_B\,\left( \Delta q_f(x_B,Q^2) +\Delta {\bar q}_f(x_B,Q^2)\right) \,,
\end{equation}
and one can write, to leading twist accuracy,
\begin{equation}
\label{eq:Jmu5}
S^\mu \Sigma(Q^2) = \frac{1}{M_N}\sum_f \langle P,S| {\bar \Psi}_f\gamma^\mu \gamma_5 \Psi_f |P,S\rangle \equiv \frac{1}{M_N}\langle P,S| J^\mu_5 (0) |P,S\rangle \,,
\end{equation}
where $J^\mu_5$ is the flavor isosinglet axial vector current in QCD. $\Sigma(Q^2)$ contributes to the spin sum rule for the proton and its value was  first extracted in pioneering experiments by the European Muon Collaboration (EMC)~\cite{Ashman:1989ig,Ashman:1987hv}; best current estimates from COMPASS~\citep{Alekseev:2010ub} give $2 \,\Sigma(Q^2) = 0.32 \pm  0.03({\rm stat.}) \pm 0.03({\rm syst.})$ at $Q^2= 3$ GeV$^2$, which is in good agreement with the extraction by HERMES~\cite{Airapetian:2007mh} at $Q^2=5$ GeV$^2$ of $2\, \Sigma(Q^2) = 0.330 \pm 0.011({\rm th.}) \pm 0.025({\rm exp.}) \pm 0.028({\rm evol.})$. 
This is significantly below the ``naive" quark model expectation~\cite{Jaffe:1989jz} of $2\,\Sigma(Q^2) = 0.6\pm 0.12$, which would result from the Ellis-Jaffe sum rule equating the isosinglet and octet axial vector currents; for more detailed discussions, see \cite{Aidala:2012mv,Deur:2018roz}. This is of course the famous ``spin crisis" of the proton -- and has lead to a large body~\cite{Kuhn:2008sy} of theoretical and experimental work since. 

The role of the anomaly becomes relevant for this discussion because the isosinglet axial vector current in Eq.~(\ref{eq:Jmu5}) is not conserved, satisfying the anomaly equation
\begin{equation}
\label{eq:anomaly}
\partial^\mu J_\mu^5(x) = \frac{n_f\alpha_s}{2\pi}~{\rm Tr} \left(F_{\mu\nu}(x)\tilde{F}^{\mu\nu}(x)\right) \,,
\end{equation}
where $F_{\mu\nu}$ is the QCD field strength tensor, its dual ${\tilde F}^{\mu\nu} = \frac{1}{2}\epsilon^{\mu\nu\rho\sigma} F_{\rho\sigma}$,  $n_f$ is the number of light quark flavors and $\alpha_s=\frac{g^2}{4\pi}$, where $g$ is the QCD coupling. 
One may however rewrite Eq.~(\ref{eq:Jmu5}) in terms of a conserved current as 
\begin{equation}
S^\mu \Sigma(Q^2) = S^\mu {\tilde \Sigma}(Q^2) + 2\, n_f \frac{1}{M_N}\,\langle P,S| K^\mu|P,S\rangle \,,
\label{eq:real-Sigma}
\end{equation}
where $S^\mu {\tilde \Sigma}(Q^2) = \frac{1}{M_N}\langle P,S| {\tilde J}_\mu^5 |P,S\rangle$. Here ${\tilde J}_\mu^5 = J_\mu^5 - 2 n_f K_\mu$ is a conserved current since the anomaly satisfies the equation 
\begin{equation}
\partial_\mu J^\mu_5 = 2\, n_f \,\partial_\mu K^\mu \,,
\end{equation}
with the Chern-Simons current $K^\mu$ defined to be 
\begin{equation}
\label{eq:C-S}
K_\mu = \frac{\alpha_S}{4 \pi} \epsilon_{\mu\nu\rho\sigma} \left[  A_a^\nu \left( \partial^\rho A_a^\sigma -\frac{1}{3} g f_{abc} A_b^\rho A_c^\sigma \right)  \right]\,.
\end{equation}

One possible explanation for the small value of $\Delta \Sigma$, advanced early \cite{Altarelli:1988nr,Carlitz:1988ab,Lampe:1998eu} after the EMC discovery, is that if first moment of $g_1$ were providing information on ${\tilde \Sigma}$ rather than $\Sigma$, that would provide a potential resolution of the spin crisis with the framework of the parton model itself. More specifically, it was argued on the basis of the gauge structure of $K^\mu$ that one could write 
\begin{equation}
 {\tilde \Sigma}(Q^2) = \Sigma(Q^2) - \frac{n_f\alpha_S}{2\pi}  \Delta G \,,
\end{equation}
where $\Delta G$ is the gluon helicity pdf. If $\Delta G$ is large, this would provide a natural explanation of the spin crisis. However as pointed out by Jaffe and Manohar, this identification is 
intrinsically problematic because while the gluon helicity pdf $\Delta G$ is manifestly gauge invariant, the same cannot be said of the Chern-Simons current. 
The latter is not gauge invariant under large gauge transformations $U$, which give,
\begin{equation}
K_\mu \rightarrow K_\mu + i\frac{g}{8 \pi^2}\epsilon_{\mu\nu\alpha\beta}\partial^\nu \left((U^\dagger \partial^\alpha U) A^\beta\right) + \frac{1}{24\pi^2} \epsilon_{\mu\nu\alpha\beta}\left[(U^\dagger \partial^\nu U) (U^\dagger \partial^\alpha U) (U^\dagger \partial^\beta U)\right]\,.
\end{equation}

The resolution of the problem, as discussed by Jaffe and Manohar~\cite{Jaffe:1989jz} (see also \cite{Efremov:1989sn}), lies in how one takes limits when $U_A(1)$ is broken. This is because the breaking of this symmetry lifts an apparent pole in the forward scattering amplitude. Indeed this is the fundamental reason why the $\eta^\prime$-meson gets a mass (distinct from the pseudoscalar octet) in QCD~\cite{Witten:1979vv,Veneziano:1979ec,Bass:2004xa,Shore:2007yn}. We will now spell out the argument as sketched in \cite{Jaffe:1989jz}. For our convenience, and that of the reader familiar with their paper, we will use their notations. 

We begin by first considering the r.h.s of Eq.~(\ref{eq:Jmu5}) and writing its off-forward counterpart as
\begin{equation}
\label{eq:off-forward-Jmu5}
\frac{1}{M_N}\langle P^\prime,S| J^\mu_5 (0)|P,S\rangle = \Sigma(Q^2,t)\,S^\mu + h(Q^2,t)\,l\cdot S \,l^\mu \,,
\end{equation}
where $l^\mu = {P^\prime}^\mu - P^\mu$ is the momentum transfer between the outgoing and incoming proton and $t=l^2$. Here $\Sigma(Q^2,t)$ and $h(Q^2,t)$ can be interpreted as form factors that represent respectively the coupling of the isosinglet axial vector charge and the isosinglet pseudoscalar charge to the proton at finite momentum transfer. In particular, the former represents the triangle diagram of the anomaly 
\begin{equation}
i l\cdot S\, {\kappa}(Q^2,t) = \frac{1}{M_N}\langle P^\prime,S| \frac{\alpha_s n_f}{2\pi}~{\rm Tr} \left(F_{\mu\nu}\tilde{F}^{\mu\nu}\right)(0)|P,S\rangle \,,
\end{equation}
which, as suggested by the r.h.s, represents the coupling of topological charge to the nucleon. The other form factor $h(t)$ represents the isosinglet pseudoscalar form factor, given by the coupling of the $\eta^\prime$ meson to the nucleon~\cite{Veneziano:1989ei,Hatsuda:1989bi,Efremov:1989sn}. 

Then from the anomaly equation (Eq.~(\ref{eq:anomaly})), and from Eq.~(\ref{eq:off-forward-Jmu5}), one obtains,
\begin{equation}
\kappa(t) = \Sigma(Q^2,t) + t \,h(Q^2,t) \, .
\end{equation}
Further, since the $\eta^\prime$ is massive, $h(t)$ has no pole, which gives $\Sigma(Q^2,0)=\kappa(Q^2,0)$. However, as stated in \cite{Jaffe:1989jz}, the triangle graph gives $S^\mu\,\Sigma(Q^2,t) \propto -i \frac{\alpha_s}{2\pi}\, \frac{l^\mu}{l^2}\, {\rm  Tr}\,(F{\tilde F})$. One cannot therefore naively take the forward limit. 

More specifically, the statement that ``$h(t)$ has no pole" and the observation that the triangle graph has an infrared pole are intimately connected and it is the interplay between the two that leads to a finite result. 
Indeed, as noted in \cite{Jaffe:1989jz}, the limit of momentum zero transfer must then be understood by writing the r.h.s of Eq.~(\ref{eq:off-forward-Jmu5}) as 
\begin{equation}
S^\mu\, \Sigma(Q^2,t) + l\cdot S\, l^\mu\, h(Q^2,t) \longrightarrow \frac{l\cdot S\,l^\mu}{t} {\kappa}(Q^2,t) + \left(S^\mu -\frac{l\cdot S\, l^\mu}{t}\right) \lambda(Q^2,t) \,.
\end{equation}
This decomposition separates the triangle graph from other contributions. For the forward matrix element of $J_\mu^5$ to appear as a smooth limit of $l^\mu\rightarrow 0$ (because the pole must be lifted by the mass of the $\eta^\prime$)  the 
triangle contribution  must cancel  a similar contribution extracted from the pseudoscalar coupling, as suggested by the second term above. One way to think about this is that there is a mixing of the contributions from the topological charge and an isosinglet component of a pseudo-Goldstone nonet which can be separated out in this manner~\cite{Diakonov:1981nv,Shore:1990zu,Ji:1990fj,Shore:1991dv,Efremov:1990nj,Dvali:2005an,Dvali:2005ws,Kharzeev:2015xsa}. For this to hold, one must of course 
have $\lambda(0) = \kappa(0)$.  As we will discuss in paper II, this follows from the Wess-Zumino-Witten 
term~\cite{Wess:1971yu,Witten:1983tw,Leutwyler:1996sa} for the $\eta^\prime$. 

What survives then on the r.h.s as $l^\mu\rightarrow 0$ is  $S^\mu \lambda(0)\equiv S^\mu \kappa(0)$. This gives the result
\begin{equation}
 \Sigma(Q^2) = \frac{n_f\,\alpha_s}{2\pi\,M_N} 
 \lim_{l_\mu \to 0}\langle P^\prime,S|\frac{ 1 }{i l\cdot s  } {\rm Tr} \left(F\tilde{F}\right)(0) |P,S\rangle\,\,.
\label{eq:sigma-anomaly}
\end{equation}
While the matrix element of ${\rm Tr} \left(F\tilde{F}\right)(0)$ is naively zero in the forward limit, the matrix element as defined above is finite when one combines the contribution from the density matrix $|P^\prime\rangle\langle P|$ and the triangle operator. This is often done in careful perturbative QCD computations by introducing a mass term or like infrared regulator which cancels between the two to give the finite result~\cite{Manohar:1990jx,Bodwin:1990fk,Vogelsang:1990ug}.  However they do not arrive at the expression in Eq.~(\ref{eq:sigma-anomaly}) because they do not further impose the constraints required by the soft dynamics\footnote{For early discussions of these in the context of $U_A(1)$ breaking by instantons, following the seminal work of t'Hooft~\cite{tHooft:1976snw,tHooft:1986ooh}, we refer the reader to Refs.~\cite{Forte:1989qq,Forte:1990xb,Dorokhov:1993ym,Qian:2015wyq}. We note that while instantons provide an attractive dynamical mechanism, this interpretation of the phenomenon is by no means unique.} of $U_A(1)$ breaking in QCD. If they are not imposed, anomaly matching cannot occur and there  will remain an unrequited pole from a pseudoscalar coupling to the triangle graph~\cite{Preskill:1990fr}. For elegant reviews of the role of the anomaly and the  Wess-Zumino-Witten term for the $\eta^\prime$ in chiral effective lagrangians for the pseudoscalar nonet, we refer the reader to \cite{Leutwyler:1997yr,HerreraSiklody:1996pm}. Lattice computations of $\Sigma(Q^2)$ implementing anomalous Ward identities are discussed in \cite{Liang:2018pis} and references therein. 

As we noted previously, our purpose here is to go beyond the discussion of the triangle anomaly in $\Sigma(Q^2)$ {\it a la} Eq.~(\ref{eq:Jmu5}) and to discuss its role in $g_1(x,Q^2)$ itself in both Bjorken and Regge asymptotics. We obtain the striking result that the formal structure of our results is identical in the two asymptotic limits of perturbative QCD. Our quantitative results for the latter in particular are novel.  As the discussion above suggests,  they have strong implications for our understanding spin diffusion at small $x_B$; these  
 will be discussed at length in Papers II \& III. 

The outline of the paper is as follows. In section~\ref{sec:world-line-hadron-tensor}, we will extend the worldline formalism developed for unpolarized DIS by us previously~\cite{Tarasov:2019rfp} to the case of polarized DIS. We will first write down the most general expression for the 
box diagram corresponding to $g_1(x_B,Q^2)$. We will then consider the Bjorken asymptotics of $Q^2\rightarrow\infty$ in section~\ref{sec:Bj} and demonstrate explicitly in Sec.~\ref{subsec:Bj-worldline} how the triangle anomaly appears and  provides the leading contribution in this asymptotics. Our result is given in Eq.~(\ref{eq:g1-Bj-final}). A discussion of this result and its implications is provided in Sec.~\ref{subsec:Bj-discussion} for readers who may not be interested in the details of the worldline derivation. Our result here is of interest firslyt because most treatments in the literature are of $\Sigma(Q^2)$ rather than $g_1(x_B,Q^2)$ itself. Moreover, unlike these discussions, we do not make use of the OPE. Our worldline framework allows us to classify graphs into those that contain the anomaly structure, and those that do not, with the latter being suppressed in the Bjorken limit. Such a  classification may be of value in the computations of other DIS observables. We also comment on the consistency of our results with an analysis of the perturbative evolution of $\Sigma (Q^2)$ and $\Delta G(Q^2)$ to high loop accuracy. 

In section~\ref{sec:Regge}, we show that the anomaly provides the leading contribution in the Regge asymptotics of $x_B\rightarrow 0$. The worldline derivation in Sec.~\ref{subsec:Rj-worldline} is very similar to that of the previous section.  Though our final result (Eq.~(\ref{resultRj})) is formally identical to that in the Bjorken limit, subtle differences in the two derivations indicate that the computation of the matrix element of the anomaly will differ both  qualitatively and quantitatively in the two limits. This is discussed in Sec.~\ref{subsec:Rj-discussion}.

A final section summarizes our results and briefly discusses their implications for the computation of $g_1(x_B,Q^2)$. As noted earlier, this computation in the Regge limit will be discussed at length in Papers II\,\& III. Appendix A discusses details of the computation of the box diagram of polarized DIS in the worldline formalism. Appendix B provides detailed expressions for coefficient functions encountered in intermediate steps of the computation. The  computation of the triangle graph in this formalism is discussed in Appendix C.  


\section{Worldline representation of antisymmetric part of the hadron tensor in polarized DIS}
\label{sec:world-line-hadron-tensor}
To compute $g_1(x,Q^2)$, we will require (see Eq.~(\ref{WA})) the antisymmetric part $\tilde{W}^{\mu\nu}$ of the hadron tensor in the worldline representation of DIS introduced in \cite{Tarasov:2019rfp}. One first reexpresses the hadron tensor in Eq.~(\ref{hT}) in terms of the second derivative of the effective action $\Gamma[A]$ with respect to the electromagnetic field $a_\mu(x)$  corresponding to the insertion of  incoming virtual photons $\gamma^\ast$ at two spacetime points: 
\begin{eqnarray}
W^{\mu\nu}(q, P, S) = \frac{1}{\pi e^2}{\rm Im}\  \int d^4x~ e^{iqx}\langle P,S| \frac{\delta^2 \Gamma[a, A]}{\delta a_\mu(\frac{x}{2})\delta a_\nu(-\frac{x}{2})} |P,S\rangle\,.
\label{TprodEffex}
\end{eqnarray}
Here $A$ denotes the gluon background field of the target. 

It is sufficient for our discussion of the triangle graph to work with the one loop QED+QCD representation\footnote{For a discussion of higher loop contributions to the effective action, we refer the reader to \cite{Fliegner:1997ra,Pawlowski:2008xh,Magnea:2013lna}. } of the worldline  effective action \cite{Schubert:2001he},
\begin{eqnarray}
&&\Gamma_{QCD}[a, A] = -\frac{1}{2} \int^T_0 \frac{dT}{T} {\rm Tr_c} \int \mathcal{D}x \int \mathcal{D} \psi \exp\Big\{-\int^T_0 d\tau \Big(\frac{1}{4} \dot{x}^2 + \frac{1}{2}\psi_\mu\dot{\psi}^\mu + ig\dot{x}^\mu (A_\mu + a_\mu) - ig \psi^\mu \psi^\nu F_{\mu\nu}(A+a)\Big)\Big\}\,,
\nonumber\\
\label{MLag}
\end{eqnarray}
which is characterized by 0+1 dimensional worldline trajectories of  Boson ($x^\mu(\tau)$) and Grassmann ($\psi^\mu(\tau)$) variables coupled to background electromagnetic ($a_\mu$) and gluon ($A_\mu$) fields. Note that the boson functional integral has periodic (P) boundary conditions while the Grassmann functional integral has anti-periodic (AP) boundary conditions. 

It is convenient to rewrite the operator in Eq. (\ref{TprodEffex}) in terms of the Fourier transformation of the effective action,
\begin{eqnarray}
{\tilde{\Gamma}}^{\mu\nu}[k_1, k_3] \equiv \int d^4z_1 d^4z_3 \frac{\delta^2 \Gamma[a, A]}{\delta a_\mu(z_1)\delta a_\nu(z_3)}|_{a=0} \,e^{ik_1 z_1} e^{ik_3 z_3}\,,
\label{secder}
\end{eqnarray}
where $k_1$ and $k_3$ denote the incoming photon four-momenta; separating out the antisymmetric part of Eq. (\ref{TprodEffex}), we obtain\footnote{Note that Eq. (\ref{WModFu}) is written in Euclidean space-time with signature $\eta = (1,1,1,1)$. In our calculation, we perform an analytical continuation to Minkowski space-time with $g = (1, -1, -1, -1)$ by the replacement $\eta_{\mu\nu}\to - g_{\mu\nu}$, see Ref. \cite{Schubert:2001he}.}
\begin{eqnarray}
i\tilde{W}^{\mu\nu}(q, P, S) = \frac{ 1 }{ 2 \pi e^2}\,{\rm Im}\ \int d^4x \,e^{-iqx} 
 \int \frac{d^4k_1}{(2\pi)^4} \int \frac{d^4k_3}{(2\pi)^4} e^{-ik_1 \frac{x}{2}} e^{ik_3 \frac{x}{2} }\langle P,S| {\tilde{\Gamma}}^{\mu\nu}_A[k_1, k_3] |P,S\rangle\,,
\label{WModFu}
\end{eqnarray}
where ${\tilde{\Gamma}}^{\mu\nu}_A[k_1, k_3] \equiv {\tilde{\Gamma}}^{\mu\nu}[k_1, k_3] - (\mu\leftrightarrow\nu)$. 

The hadron tensor in Eq. (\ref{TprodEffex}) is taken in the forward limit when $k_1=-k_3=-q$. However as discussed in the introduction, to obtain the infrared pole of the anomaly, one needs to calculate the off-forward matrix element $\langle P'|\dots|P\rangle$ in Eq. (\ref{TprodEffex}), where $P'-P\equiv l$, and then take the limit $l\to0$ in the final expression. Hence the incoming photon momenta in our computation of $\Gamma^{\mu\nu}_A[k_1, k_3]$ are kept distinct.

Taking the second derivative of the effective action, we obtain
\begin{eqnarray}
{\tilde{\Gamma}}^{\mu\nu}_A[k_1, k_3] &=& \frac{e^2 e_f^2 }{2}\int^\infty_0\frac{dT}{T} ~{\rm Tr_c} \int \mathcal{D}x\int \mathcal{D}\psi 
\Big[ V^\mu_1(k_1)V^\nu_3(k_3) - (\mu\leftrightarrow\nu)\Big] \label{photamp}\\
&\times& \exp\Big\{-\int^T_0 d\tau \Big(\frac{1}{4} \dot{x}^2 + \frac{1}{2}\psi_\mu\dot{\psi}^\mu + ig\dot{x}^\mu A_\mu - ig \psi^\mu \psi^\nu F_{\mu\nu}\Big)\Big\}\,,
\nonumber
\end{eqnarray}
where
\begin{eqnarray}
V^\mu_i(k_i) \equiv \int^T_0 d\tau_i  (\dot{x}^\mu_i + 2i \psi^\mu_i k_j\cdot\psi_j ) e^{ik_i\cdot x_i}\,,
\end{eqnarray}
is the vertex corresponding to the interaction of a  worldline  with the external electromagnetic current, and $x_i\equiv x(\tau_i)$, $\psi_i\equiv \psi(\tau_i)$.

Expanding the worldline action up to the second order in the background field\footnote{In discussions of the anomaly, it is often convenient to impose Fock-Schwinger gauge $x\cdot A=0$, and expand the result in powers of $F_{\mu\nu}$-- see~\cite{Mueller:2017arw} for instance, for an explicit derivation of the anomaly equation in the worldline formalism. For our discussion of the anomaly, the approach here is sufficient; that it is so is a non-trivial feature of the non-Abelian axial anomaly~\cite{AlvarezGaume:2005qb,Bilal:2008qx}.}, one obtains, 
\begin{eqnarray}
&&{\tilde{\Gamma}}^{\mu\nu}_A[k_1, k_3] = (-ig)^2 \frac{e^2 e_f^2 }{2}\int^\infty_0\frac{dT}{T} ~{\rm Tr_c} \int \mathcal{D}x\int \mathcal{D}\psi~\exp\Big\{-\int^T_0 d\tau \Big(\frac{1}{4} \dot{x}^2 + \frac{1}{2}\psi \cdot\dot{\psi} \Big)\Big\}
\\
&&\times~\Big[ V^\mu_1(k_1)V^\nu_3(k_3) \int^T_0 d\tau_2 \Big( \dot{x}^\alpha_2 A_\alpha(x_2) +  2 \psi^\alpha_2 \psi^\lambda_2 \partial_\lambda A_{\alpha}(x_2)\Big) \int^T_0 d\tau_4 \Big( \dot{x}^\beta_4 A_\beta(x_4) + 2 \psi^\beta_4 \psi^\eta_4 \partial_\eta A_{\beta}(x_4) \Big) - (\mu\leftrightarrow\nu)\Big]\,.
\nonumber
\end{eqnarray}
One can rewrite this further in terms of the Fourier transforms of the background gauge fields
\begin{eqnarray}
A_\alpha(x_2) = \int \frac{d^4k_2}{(2\pi)^4} e^{ik_2\cdot x_2} {\tilde A}_\alpha(k_2);\ \ \ \ \ A_\beta(x_4) = \int \frac{d^4k_4}{(2\pi)^4} e^{ik_4\cdot x_4} {\tilde A}_\beta(k_4)\,,
\label{fourbackground}
\end{eqnarray}
as 
\begin{eqnarray}
&&\Gamma^{\mu\nu}_A[k_1, k_3] = \int \frac{d^4k_2}{(2\pi)^4} \int \frac{d^4k_4}{(2\pi)^4}~
\Gamma^{\mu\nu\alpha\beta}_A[k_1, k_3, k_2, k_4]~ {\rm Tr_c}({\tilde A}_\alpha(k_2) {\tilde A}_\beta(k_4))\,,
\label{Hadgamma}
\end{eqnarray}
where
\begin{eqnarray}
&&\Gamma^{\mu\nu\alpha\beta}_A[k_1, k_3, k_2, k_4] \equiv -\frac{g^2 e^2 e_f^2 }{2}\int^\infty_0\frac{dT}{T} ~ \int \mathcal{D}x\int \mathcal{D}\psi~\exp\Big\{-\int^T_0 d\tau \Big(\frac{1}{4} \dot{x}^2 + \frac{1}{2}\psi\cdot\dot{\psi} \Big)\Big\}
\nonumber\\
&&\times \Big[V^\mu_1(k_1)V^\nu_3(k_3)V^\alpha_2(k_2) V^\beta_4(k_4) - (\mu\leftrightarrow\nu)\Big]\,,
\end{eqnarray}
corresponds to the well-known box diagram of DIS with four incoming momenta $k_i$ shown in Fig. \ref{fig:boxgraph}.

\begin{figure}[htb]
 \begin{center}
 \includegraphics[width=90mm]{box}
 \end{center}
 \caption{\label{fig:boxgraph}The box diagram $\Gamma^{\mu\nu\alpha\beta}_A[k_1, k_3, k_2, k_4]$ for polarized DIS.}
 \end{figure}

Taking the product of the four worldline interaction vertices $V^\mu_i(k_i)$, and removing the terms proportional to $ \dot{x}^\mu_1\dot{x}^\nu_3$ and $\dot{x}^\alpha_2 \dot{x}^\beta_4$, which do not contribute to the antisymmetric part of the hadron tensor\footnote{Indeed, the interaction of the worldline with external particles through $\dot{x}^\mu_1\dot{x}^\nu_3$ and $\dot{x}^\alpha_2 \dot{x}^\beta_4$ coincides with that of  scalar QED. As a result, it doesn't generate any spin dependent effect.},
we obtain
\begin{eqnarray}
&&\Gamma^{\mu\nu\alpha\beta}_A[k_1, k_3, k_2, k_4] = -\frac{g^2e^2 e_f^2 }{2}\int^\infty_0\frac{dT}{T} ~ \int \mathcal{D}x\int \mathcal{D}\psi~\exp\Big\{-\int^T_0 d\tau \Big(\frac{1}{4} \dot{x}^2 + \frac{1}{2}\psi\cdot \dot{\psi} \Big)\Big\}
\nonumber\\
&&\times \prod^4_{k=1}\int^T_0 d\tau_k~ \Big[\sum^9_{n=1}\mathcal{C}^{\mu\nu\alpha\beta}_{n;(\tau_1,\tau_2,\tau_3,\tau_4)}[k_1, k_3, k_2, k_4] - (\mu\leftrightarrow\nu)\Big]e^{i\sum^4_{i=1}k_i x_i} \,.
\label{GammawithC}
\end{eqnarray}
where the coordinate ($x_i \equiv x(\tau_i)$) and Grassmann variables ($\psi_i\equiv \psi(\tau_i)$) in the coefficients $\mathcal{C}^{\mu\nu\alpha\beta}_{n;(\tau_1,\tau_2,\tau_3,\tau_4)}[k_1, k_3, k_2, k_4]$ depend on the proper time coordinates of the interaction of the worldlines with the external electromagnetic and gauge fields:
\begin{eqnarray}
&&\mathcal{C}^{\mu\nu\alpha\beta}_{1;(\tau_1,\tau_2,\tau_3,\tau_4)}[k_1, k_3, k_2, k_4] = - 4\dot{x}^\nu_3 \psi^\mu_1 \psi_1\cdot k_{1} \dot{x}^\beta_4 \psi^\alpha_2 \psi_2\cdot k_{2};
\nonumber\\
&&\mathcal{C}^{\mu\nu\alpha\beta}_{2;(\tau_1,\tau_2,\tau_3,\tau_4)}[k_1, k_3, k_2, k_4] = -4\dot{x}^\nu_3 \psi^\mu_1 \psi_1\cdot k_{1}  \dot{x}^\alpha_2 \psi^\beta_4 \psi_4\cdot k_{4};
\nonumber\\
&&\mathcal{C}^{\mu\nu\alpha\beta}_{3;(\tau_1,\tau_2,\tau_3,\tau_4)}[k_1, k_3, k_2, k_4] = -4 \dot{x}^\mu_1\psi^\nu_3 \psi_3 \cdot k_{3} \dot{x}^\alpha_2 \psi^\beta_4 \psi_4\cdot k_{4};
\nonumber\\
&&\mathcal{C}^{\mu\nu\alpha\beta}_{4;(\tau_1,\tau_2,\tau_3,\tau_4)}[k_1, k_3, k_2, k_4] = -4 \dot{x}^\mu_1\psi^\nu_3 \psi_3 \cdot k_{3}  \dot{x}^\beta_4 \psi^\alpha_2 \psi_2\cdot k_{2}
\nonumber\\
&&\mathcal{C}^{\mu\nu\alpha\beta}_{5;(\tau_1,\tau_2,\tau_3,\tau_4)}[k_1, k_3, k_2, k_4] = - 8i\dot{x}^\nu_3 \psi^\mu_1 \psi_1\cdot k_{1} \psi^\alpha_2 \psi_2\cdot k_{2} \psi^\beta_4 \psi_4\cdot k_{4};
\nonumber\\
&&\mathcal{C}^{\mu\nu\alpha\beta}_{6;(\tau_1,\tau_2,\tau_3,\tau_4)}[k_1, k_3, k_2, k_4] = - 8i \dot{x}^\mu_1\psi^\nu_3 \psi_3 \cdot k_{3} \psi^\alpha_2 \psi_2\cdot k_{2} \psi^\beta_4 \psi_4\cdot k_{4}
\nonumber\\
&&\mathcal{C}^{\mu\nu\alpha\beta}_{7;(\tau_1,\tau_2,\tau_3,\tau_4)}[k_1, k_3, k_2, k_4] = - 8i \dot{x}^\beta_4 \psi^\alpha_2 \psi_2\cdot k_{2} \psi^\mu_1 \psi_1\cdot k_{1} \psi^\nu_3 \psi_3\cdot k_{3}  ;
\nonumber\\
&&\mathcal{C}^{\mu\nu\alpha\beta}_{8;(\tau_1,\tau_2,\tau_3,\tau_4)}[k_1, k_3, k_2, k_4] = - 8i \dot{x}^\alpha_2 \psi^\beta_4 \psi_4\cdot k_{4} \psi^\mu_1 \psi_1\cdot k_{1} \psi^\nu_3 \psi_3\cdot k_{3} 
\nonumber\\
&&\mathcal{C}^{\mu\nu\alpha\beta}_{9;(\tau_1,\tau_2,\tau_3,\tau_4)}[k_1, k_3, k_2, k_4] = 16 \psi^\mu_1 \psi_1\cdot k_{1} \psi^\nu_3 \psi_3\cdot k_{3} \psi^\alpha_2 \psi_2\cdot k_{2} \psi^\beta_4 \psi_4\cdot k_{4}
\end{eqnarray}

We can further rewrite Eq.~(\ref{GammawithC}) as
\begin{eqnarray}
&&\Gamma^{\mu\nu\alpha\beta}_A[k_1, k_3, k_2, k_4]
\nonumber\\
&&= -\frac{g^2e^2 e_f^2 }{2} (2\pi)^4\delta^{(4)}(\sum^4_{i=1} k_i)\int^\infty_0\frac{dT}{T} 
 \prod^4_{k=1}\int^T_0 d\tau_k~\frac{1}{4\pi^2 T^2} \Big(\sum^9_{n=1}\mathcal{I}^{\mu\nu\alpha\beta}_{n;(\tau_1,\tau_2,\tau_3,\tau_4)}[k_1, k_3, k_2, k_4]\Big)\langle e^{i\sum^4_{i=1}k_i x_i}\rangle\,,
\label{GammawithI}
\end{eqnarray}
where
\begin{eqnarray}
&&\frac{1}{4\pi^2 T^2}~\mathcal{I}^{\mu\nu\alpha\beta}_{n;(\tau_1,\tau_2,\tau_3,\tau_4)}[k_1, k_3, k_2, k_4] \langle e^{i\sum^4_{i=1}k_i x_i}\rangle 
\nonumber\\
&&= \int \mathcal{D}x\int \mathcal{D}\psi~\Big(\mathcal{C}^{\mu\nu\alpha\beta}_{n;(\tau_1,\tau_2,\tau_3,\tau_4)}[k_1, k_3, k_2, k_4] e^{i\sum^4_{i=1}k_i x_i} - (\mu\leftrightarrow\nu)\Big)~\exp\Big\{-\int^T_0 d\tau \Big(\frac{1}{4} \dot{x}^2 + \frac{1}{2}\psi\cdot \dot{\psi} \Big)\Big\} \,.
\label{CtoI}
\end{eqnarray} 
The correlator of the exponential factors in Eq. (\ref{GammawithI}) is obtained as an intermediate step in the computation of the functional integrals in Eq. (\ref{GammawithC}); it can be expressed as 
\begin{eqnarray}
\label{eq:correlator}
&&\langle e^{i\sum^4_{i=1}k_i x_i}\rangle = \exp\Big[
k_{1} \cdot k_{2} G_B(\tau_1, \tau_2) + k_{1} \cdot k_{3} G_B(\tau_1, \tau_3) + k_{1}\cdot k_{4} G_B(\tau_1, \tau_4)
\nonumber\\
&& + k_{2} \cdot k_{3} G_B(\tau_2, \tau_3) + k_{2} \cdot k_{4} G_B(\tau_2, \tau_4) + k_{3} \cdot k_{4} G_B(\tau_3, \tau_4) \Big]\,,
\end{eqnarray}
where  
\begin{eqnarray}
G_B(\tau_i, \tau_j) = |\tau_i - \tau_j| - \frac{(\tau_i - \tau_j)^2}{T}\,,
\end{eqnarray}
is the bosonic worldline propagator~\cite{Strassler:1992zr} on a closed loop of period ${\rm T}$.
This remarkably simple result for the correlator is a generic feature of worldline path integrals and  follows from performing Wick contractions of bosonic and Grassmann worldline propagagators employing techniques pioneered by Bern and Kosower~\cite{Bern:1990ux,Bern:1991aq,Bern:1992ad} and discussed at length in \cite{Schubert:2001he}. 
The explicit expression for one of the coefficients ($\mathcal{I}^{\mu\nu\alpha\beta}_{1;(\tau_1,\tau_2,\tau_3,\tau_4)}[k_1, k_3, k_2, k_4]$) is worked out in Appendix \ref{sec:I1}. 

It is convenient to introduce a reparametrization $\tau = u T$ of the proper time variables, where $u\in [0, 1]$. With this reparametrization, Eq. (\ref{GammawithI}) can be rewritten as 
\begin{eqnarray}
\label{eq:Master}
&&\Gamma^{\mu\nu\alpha\beta}_A[k_1, k_3, k_2, k_4]
\nonumber\\
&&= -\frac{g^2e^2 e_f^2 }{8\pi^2}(2\pi)^4\delta^{(4)}(\sum^4_{i=1} k_i)
 \prod^4_{k=1}\int^1_0 du_k \Big(\sum^9_{n=1}\mathcal{I}^{\mu\nu\alpha\beta}_{n;(u_1,u_2,u_3,u_4)}[k_1, k_3, k_2, k_4]\Big)\int^\infty_0dT~ T \exp\Big[ T\Big( k_{1} \cdot k_{2} G_B(u_1, u_2)
 \nonumber\\
 &&+ k_{1} \cdot k_{3} G_B(u_1, u_3) + k_{1}\cdot k_{4} G_B(u_1, u_4) + k_{2} \cdot k_{3} G_B(u_2, u_3) + k_{2} \cdot k_{4} G_B(u_2, u_4) + k_{3} \cdot k_{4} G_B(u_3, u_4) \Big) \Big]\,.
\end{eqnarray}
The integration over the worldline period $T$ can now be performed easily, and one obtains, 
\begin{eqnarray}
&&\Gamma^{\mu\nu\alpha\beta}_A[k_1, k_3, k_2, k_4]
\nonumber\\
&&= -\frac{g^2e^2 e_f^2 }{8\pi^2}(2\pi)^4\delta^{(4)}(\sum^4_{i=1} k_i)
 \prod^4_{k=1}\int^1_0 du_k \Big(\sum^9_{n=1}\mathcal{I}^{\mu\nu\alpha\beta}_{n;(u_1,u_2,u_3,u_4)}[k_1, k_3, k_2, k_4]\Big)\Big[ - k_{1} \cdot k_{2} G_B(u_1, u_2)
 \nonumber\\
 &&- k_{1} \cdot k_{3} G_B(u_1, u_3) - k_{1}\cdot k_{4} G_B(u_1, u_4) - k_{2} \cdot k_{3} G_B(u_2, u_3) - k_{2} \cdot k_{4} G_B(u_2, u_4) - k_{3} \cdot k_{4} G_B(u_3, u_4) \Big]^{-2}\,,
 \label{GammawithIinu}
\end{eqnarray}

This result for the box diagram is noteworthy because the functional dependence on the external momenta of the gauge fields does not rely on any kinematic assumptions. Further, since the result is insensitive to color or flavor, at this level, it only depends on these external momenta and not on whether the gauge field corresponding to a momentum label is a gluon or a photon. As we will show, this property of the box diagram in the worldline formalism will prove extremely useful. In particular, in the following two sections, we will explore  the structure of the box diagram represented by Eq. (\ref{GammawithIinu}) in the physically interesting Bjorken and Regge asymptotics respectively.  We will show explicitly in both limits that the leading contribution to the box diagram is given by the triangle anomaly.

\begin{figure}[htb]
 \begin{center}
 \includegraphics[width=120mm]{two_limits}
 \end{center}
 \caption{\label{fig:two_limits}Two ``triangle" limits of the box diagram: a) the Bjorken limit, given by  Eq.~(\ref{GammawithIBj}),  b) the Regge limit given by  Eq.~(\ref{GammawithIRj}).}
 \end{figure}

In the Bjorken limit of QCD, when the virtuality of the incoming photons $Q^2\to\infty$ and $x_{B} = \frac{Q^2}{2P\cdot q}$ is fixed, the distance between the points of interaction of the worldline with the incoming photons ($\tau_1\rightarrow \tau_3$) is defined by a negligibly small number $u_1 - u_3 \sim \Lambda_{\rm QCD}^2/Q^2$, where $\Lambda_{\rm QCD}\approx 200$ MeV is the intrinsic non-perturbative scale of the theory. This limit of the box diagram is illustrated in  Fig.~\ref{fig:two_limits}a. In this limit, 
\begin{eqnarray}
&&\Gamma^{\mu\nu\alpha\beta}_A[k_1, k_3, k_2, k_4]\Big|_{Q^2\to \infty}
\nonumber\\
&&= -\frac{g^2e^2 e_f^2 }{8\pi^2}(2\pi)^4\delta^{(4)}(\sum^4_{i=1} k_i)
 \prod^4_{k=1}\int^1_0 du_k \Big(\sum^9_{n=1}\mathcal{I}^{\mu\nu\alpha\beta}_{n;(u_1,u_2,u_3,u_4)}[k_1, k_3, k_2, k_4]\Big)\Big|_{u_1 = u_3}\Big[ - k_{1} \cdot k_{2} G_B(u_1, u_2)
 \nonumber\\
 &&- k_{1} \cdot k_{3} G_B(u_1, u_3) - k_{1}\cdot k_{4} G_B(u_1, u_4) - k_{2} \cdot k_{3} G_B(u_2, u_3) - k_{2} \cdot k_{4} G_B(u_2, u_4) - k_{3} \cdot k_{4} G_B(u_3, u_4) \Big]^{-2}\,.
\label{GammawithIBj}
\end{eqnarray}
Corrections to this formula are suppressed by a relative power $1/Q^2$.

In a similar fashion, the Regge limit of perturbative QCD (pQCD) is characterized by a fixed virtuality $Q^2\gg \Lambda_{\rm QCD}^2$ and $x_{Bj}\to 0$. In these asymptotics, the interaction of the worldline with the background gluons corresponds to an instantaneous interaction with a shock wave; in the box diagram, this corresponds to $\tau_2\rightarrow \tau_4$, or equivalently, $u_2\simeq u_4$, as shown in Fig.~\ref{fig:two_limits}b. As a result, in this limit,
\begin{eqnarray}
&&\Gamma^{\mu\nu\alpha\beta}_A[k_1, k_3, k_2, k_4]\Big|_{x_{Bj}\to 0}
\nonumber\\
&&= -\frac{g^2e^2 e_f^2 }{8\pi^2}(2\pi)^4\delta^{(4)}(\sum^4_{i=1} k_i)
 \prod^4_{k=1}\int^1_0 du_k \Big(\sum^9_{n=1}\mathcal{I}^{\mu\nu\alpha\beta}_{n;(u_1,u_2,u_3,u_4)}[k_1, k_3, k_2, k_4]\Big)\Big|_{u_2 = u_4}\Big[ - k_{1} \cdot k_{2} G_B(u_1, u_2)
 \nonumber\\
 &&- k_{1} \cdot k_{3} G_B(u_1, u_3) - k_{1}\cdot k_{4} G_B(u_1, u_4) - k_{2} \cdot k_{3} G_B(u_2, u_3) - k_{2} \cdot k_{4} G_B(u_2, u_4) - k_{3} \cdot k_{4} G_B(u_3, u_4) \Big]^{-2}\,.
\label{GammawithIRj}
\end{eqnarray}
In analogy to the Bjorken limit, corrections to this expression are suppressed by a relative power $\sim Q_s^2/M^2$, where $M^2= 2\,xP\cdot q$ is a large scale when $s\rightarrow \infty$. Here $x$ denotes the longitudinal momentum fraction of the proton momentum carried by the background gluon, and $Q_s$ denotes the typical transverse momentum\footnote{We will discuss this emergent scale further in section~\ref{sec:Regge}.} of the background gluons in Regge asymptotics. 

As we noted earlier, the box diagram is only sensitive to the external momentum labels of the external gauge fields. It  will therefore have an identical structure when any two of the proper time values that they correspond to (see Eq.~(\ref{eq:correlator})) are set equal to each other.  Hence the expressions for the Bjorken limit (Eq.~(\ref{GammawithIBj})) and the Regge limit (Eq.~(\ref{GammawithIRj})) can both be understood simply in this formalism as proper time limits of the box diagram;  indeed, as sketched in Fig.~\ref{fig:two_limits}, the triangle structure has a clear visual representation in both asymptotics. This is not the case for the usual OPE language of pQCD where the former is manifest while the latter is not. 

An interesting consequence of our results will be that the structure of the matrix elements for $g_1(x_{\rm Bj},Q^2)$ will be identical in both Bjorken and Regge asymptotics. However the underlying physics of the matrix elements, in particular their QCD evolution, will be quite different in the two kinematics limits. 


\section{The triangle anomaly in the Bjorken limit of the box diagram}
\label{sec:Bj}
In this section, we will compute the box diagram of Fig. \ref{fig:boxgraph} in the Bjorken limit. This corresponds to a resolution scale corresponding to a transverse area in the 
proton which vanishes with $1/Q^2\rightarrow \infty$; in the QCD worldline formalism, this corresponds to $u_1 \simeq u_3$ in Eq.~(\ref{eq:Master}), resulting in Eq.~(\ref{GammawithIBj}). 
\subsection{Worldline computation of box diagram in the Bjorken limit}
\label{subsec:Bj-worldline}
 We have explicitly computed the first of the coefficients $\mathcal{I}^{\mu\nu\alpha\beta}_{n;(u_1,u_2,u_3,u_4)}[k_1, k_3, k_2, k_4]$ explicitly in Appendix A; the rest can be computed similarly. The expressions for all of the nine terms are provided in Appendix B. The sum of these, for $u_1=u_3$, is given by  
\begin{eqnarray}
&&\Big(\sum^9_{n=1}\mathcal{I}^{\mu\nu\alpha\beta}_{n;(u_1,u_2,u_3,u_4)}[k_1, k_3, k_2, k_4]\Big)\Big|_{u_1 = u_3} = \epsilon^{\mu\nu\eta}_{\ \ \ \ \kappa} (k_{1\eta} - k_{3\eta})
\nonumber\\
&&\times \Big( k_2 \cdot k_4 \Big[ \dot{G}^2_B(u_1, u_4) - \dot{G}_B(u_1, u_4) \dot{G}_B(u_2, u_4) + \dot{G}_B(u_1, u_2) \big(\dot{G}_B(u_1, u_2) + \dot{G}_B(u_2, u_4) + \dot{G}_B(u_4, u_1) \big)  \Big] \epsilon^{\kappa\alpha\beta\sigma} k_{2\sigma}
\nonumber\\
&&+ k_2\cdot k_4 \Big[ -\dot{G}^2_B(u_1, u_2) - \dot{G}_B(u_1, u_2) \dot{G}_B(u_2, u_4) + \dot{G}_B(u_1, u_4) \big(\dot{G}_B(u_1, u_2) + \dot{G}_B(u_2, u_4) + \dot{G}_B(u_4, u_1)\big) \Big] \epsilon^{\kappa\alpha\beta\sigma} k_{4\sigma}
\nonumber\\
&&+ \Big[ - \dot{G}_B(u_1, u_4) \dot{G}_B(u_2, u_4) +  \dot{G}^2_B(u_1, u_4)
+ \dot{G}_B(u_1, u_2) \big(\dot{G}_B(u_1, u_2) + \dot{G}_B(u_2, u_4) + \dot{G}_B(u_4, u_1)\big) \Big] \epsilon^{\kappa\alpha\sigma\lambda} k^\beta_{2} k_{2\sigma} k_{4\lambda}
\nonumber\\
&&+ \dot{G}^2_B(u_1, u_4) \epsilon^{\kappa\alpha\sigma\lambda} k^\beta_{4} k_{2\sigma} k_{4\lambda}
- \dot{G}^2_B(u_1, u_2) \epsilon^{\kappa\beta\sigma\lambda} k^\alpha_{2} k_{2\sigma} k_{4\lambda} 
\nonumber\\
&&+ \Big[ - \dot{G}_B(u_1, u_2) \dot{G}_B(u_2, u_4) - \dot{G}^2_B(u_1, u_2)
+ \dot{G}_B(u_1, u_4) \big(\dot{G}_B(u_1, u_2) + \dot{G}_B(u_2, u_4) + \dot{G}_B(u_4, u_1)\big) \Big] \epsilon^{\kappa\beta\sigma\lambda} k^\alpha_{4} k_{2\sigma} k_{4\lambda} 
\nonumber\\
&&- \epsilon^{\alpha\lambda\beta\sigma}  k^\kappa_{2} k_{2\lambda} k_{4\sigma} - \epsilon^{\alpha\lambda\beta\sigma} k^\kappa_{4} k_{2\lambda} k_{4\sigma}\Big)\,.
\label{suminBj}
\end{eqnarray}
This expression looks formidable with a large number of Lorentz structures; such structures are also obtained in the well-known perturbative computation of the triangle diagram of the anomaly~\cite{Rosenberg:1962pp,Armillis:2009sm}. We will show explicitly that this expression can be greatly simplified using a few identities and a mass-shell constraint on the background gluons.

We first use the identity\footnote{\label{cyclic} This result can be easily obtained from the standard cyclic identity $v^\mu \epsilon^{\nu\alpha\beta\sigma} + v^\nu \epsilon^{\alpha\beta\sigma\mu} + v^\alpha \epsilon^{\beta\sigma\mu\nu} + v^\beta \epsilon^{\sigma\mu\nu\alpha} + v^\sigma \epsilon^{\mu\nu\alpha\beta} = 0$, where $v^\mu$ is an arbitrary four-vector and a mass-shell condition $k^2_2 = k^2_4= 0$ is imposed on the background gluons.}
\begin{eqnarray}
\label{Levi-Civita}
\epsilon^{\alpha\lambda\beta\sigma}  k^\kappa_{2} k_{2\lambda} k_{4\sigma} = - k^\alpha_{2} \epsilon^{\lambda\beta\sigma\kappa} k_{2\lambda} k_{4\sigma}
- k^\beta_{2} \epsilon^{\sigma\kappa\alpha\lambda}  k_{2\lambda} k_{4\sigma}
- k_{2}\cdot k_{4} \epsilon^{\kappa\alpha\lambda\beta}  k_{2\lambda} \,,
\end{eqnarray}
for the last but one term in Eq. (\ref{suminBj}), and a similar identity for the last term, to reexpress it as 
\begin{eqnarray}
&&\Big(\sum^9_{n=1}\mathcal{I}^{\mu\nu\alpha\beta}_{n;(u_1,u_2,u_3,u_4)}[k_1, k_3, k_2, k_4]\Big)\Big|_{u_1 = u_3} = \epsilon^{\mu\nu\eta}_{\ \ \ \ \kappa} (k_{1\eta} - k_{3\eta})
\nonumber\\
&&\times \Big( \Big[ -1 + \dot{G}^2_B(u_1, u_4) - \dot{G}_B(u_1, u_4) \dot{G}_B(u_2, u_4) + \dot{G}_B(u_1, u_2) \big(\dot{G}_B(u_1, u_2) + \dot{G}_B(u_2, u_4) + \dot{G}_B(u_4, u_1) \big)  \Big]
\nonumber\\
&&\times( k_2 \cdot k_4 \epsilon^{\kappa\alpha\beta\sigma} k_{2\sigma} + \epsilon^{\kappa\alpha\sigma\lambda} k^\beta_{2} k_{2\sigma} k_{4\lambda})
+ \Big[ 1 - \dot{G}^2_B(u_1, u_2) - \dot{G}_B(u_1, u_2) \dot{G}_B(u_2, u_4)
\nonumber\\
&&+ \dot{G}_B(u_1, u_4) \big(\dot{G}_B(u_1, u_2) + \dot{G}_B(u_2, u_4) + \dot{G}_B(u_4, u_1)\big) \Big]
 (k_2\cdot k_4 \epsilon^{\kappa\alpha\beta\sigma} k_{4\sigma} + \epsilon^{\kappa\beta\sigma\lambda} k^\alpha_{4} k_{2\sigma} k_{4\lambda})
\nonumber\\
&&+ \Big[ -1 + \dot{G}^2_B(u_1, u_4) \Big] \epsilon^{\kappa\alpha\sigma\lambda} k^\beta_{4} k_{2\sigma} k_{4\lambda}
+ \Big[ 1 - \dot{G}^2_B(u_1, u_2) \Big] \epsilon^{\kappa\beta\sigma\lambda} k^\alpha_{2} k_{2\sigma} k_{4\lambda} \Big)\, .
\end{eqnarray}

Now recall that the $\alpha,\beta$ open indices on the l.h.s are contracted with the gauge fields, as shown in 
Eq.~(\ref{Hadgamma}).  Since the background gluons are on mass-shell, $k^\beta_{4} A_\beta(k_4) = 0$ and $k^\alpha_{2} A_\alpha(k_2) = 0$. Hence the last two tensorial structures in the above equation don't contribute. Further, the first two structures can be simplified using 
\begin{eqnarray}
k_2 \cdot k_4 \epsilon^{\kappa\alpha\beta\sigma} k_{2\sigma} + \epsilon^{\kappa\alpha\sigma\lambda} k^\beta_{2} k_{2\sigma} k_{4\lambda} \rightarrow  - k^\kappa_{2} \epsilon^{\alpha\beta\sigma\lambda}  k_{2\sigma} k_{4\lambda}\,,
\end{eqnarray}
where we have used $k^\alpha_{2} A_\alpha(k_2) = 0$ to eliminate a term in Eq.~(\ref{Levi-Civita}), and likewise, 
\begin{eqnarray}
k_2\cdot k_4 \epsilon^{\kappa\alpha\beta\sigma} k_{4\sigma} + \epsilon^{\kappa\beta\sigma\lambda} k^\alpha_{4} k_{2\sigma} k_{4\lambda} \rightarrow  k^\kappa_{4} \epsilon^{\alpha\beta\sigma\lambda}  k_{2\sigma} k_{4\lambda}\,.
\end{eqnarray}
As a result, the sum of coefficients simplifies to
\begin{eqnarray}
&&\Big(\sum^9_{n=1}\mathcal{I}^{\mu\nu\alpha\beta}_{n;(u_1,u_2,u_3,u_4)}[k_1, k_3, k_2, k_4]\Big)\Big|_{u_1 = u_3} = - \epsilon^{\mu\nu\eta}_{\ \ \ \ \kappa} (k_{1\eta} - k_{3\eta})
\\
&&\times \Big( \Big[ -1 + \dot{G}^2_B(u_1, u_4) - \dot{G}_B(u_1, u_4) \dot{G}_B(u_2, u_4) + \dot{G}_B(u_1, u_2)\, {\cal X}(u_1,u_2,u_4)   \Big] k^\kappa_{2} \epsilon^{\alpha\beta\sigma\lambda}  k_{2\sigma} k_{4\lambda}
\nonumber\\
&&+ \Big[ -1 + \dot{G}^2_B(u_1, u_2) + \dot{G}_B(u_1, u_2) \dot{G}_B(u_2, u_4) - \dot{G}_B(u_1, u_4)\, {\cal X}(u_1,u_2,u_4) \Big] k^\kappa_{4} \epsilon^{\alpha\beta\sigma\lambda}  k_{2\sigma} k_{4\lambda} \Big)\,,
\nonumber
\end{eqnarray}
where 
\begin{eqnarray}
{\cal X}(u_1,u_2,u_4) =\dot{G}_B(u_1, u_2) + \dot{G}_B(u_2, u_4) + \dot{G}_B(u_4, u_1) \,.
\end{eqnarray}
This expression can be written more compactly as 
\begin{eqnarray}
\Big(\sum^9_{n=1}\mathcal{I}^{\mu\nu\alpha\beta}_{n;(u_1,u_2,u_3,u_4)}[k_1, k_3, k_2, k_4]\Big)\Big|_{u_1 = u_3} &=& - \frac{1}{2}\epsilon^{\mu\nu\eta}_{\ \ \ \ \kappa} (k_{1\eta} - k_{3\eta})
\,\Big[ - 2 + {\cal X}^2(u_1,u_2,u_4)
\nonumber \\
&+&  \dot{G}^2_B(u_1, u_2) - \dot{G}^2_B(u_2, u_4) + \dot{G}^2_B(u_4, u_1)
 \Big] ( k^\kappa_{2} + k^\kappa_{4} ) \epsilon^{\alpha\beta\sigma\lambda}  k_{2\sigma} k_{4\lambda}\,\,.
\end{eqnarray}

One can show that 
\begin{eqnarray}
{\cal X}(u_1,u_2,u_4) \equiv \dot{G}_B(u_1, u_2) + \dot{G}_B(u_2, u_4) + \dot{G}_B(u_4, u_1) = -G_F(u_1, u_2)G_F(u_2, u_4)G_F(u_4, u_1)
\label{gfprod}
\end{eqnarray}
where $G_F(u_i,u_j) = {\rm sign}(u_i - u_j)$ is the fermionic worldline propagator. Using then ${\cal X}^2 =1$ and
\begin{eqnarray}
1 - \dot{G}^2_B(u_i, u_j) = 4 \,G_B(u_i, u_j)\,,
\label{b2dot}
\end{eqnarray}
we obtain 
\begin{eqnarray}
\Big(\sum^9_{n=1}\mathcal{I}^{\mu\nu\alpha\beta}_{n;(u_1,u_2,u_3,u_4)}[k_1, k_3, k_2, k_4]\Big)\Big|_{u_1 = u_3} &=& - 2\epsilon^{\mu\nu\eta}_{\ \ \ \ \kappa} (k_{1\eta} - k_{3\eta})
 \Big[ - G_B(u_1, u_2) + G_B(u_2, u_4) - G_B(u_4, u_1)
 \Big] \nonumber \\
 &\times& ( k^\kappa_{2} + k^\kappa_{4} ) \epsilon^{\alpha\beta\sigma\lambda}  k_{2\sigma} k_{4\lambda}\,\,.
 \label{sumcoefBj}
\end{eqnarray}

Substituting this back to Eq. (\ref{GammawithIBj}),  we obtain the following result for the box diagram in the Bjorken limit of QCD:
\begin{eqnarray}
&&\Gamma^{\mu\nu\alpha\beta}_A[k_1, k_3, k_2, k_4]\Big|_{Q^2\to \infty}
= \frac{g^2e^2 e_f^2 }{4\pi^2}\epsilon^{\mu\nu\eta}_{\ \ \ \ \kappa} (k_{1\eta} - k_{3\eta})(2\pi)^4\delta^{(4)}(\sum^4_{i=1} k_i)
 \prod^4_{k=1}\int^1_0 du_k 
 \nonumber\\
&&\times \Big[ - G_B(u_1, u_2) + G_B(u_2, u_4) - G_B(u_4, u_1)
 \Big] ( k^\kappa_{2} + k^\kappa_{4} ) \epsilon^{\alpha\beta\sigma\lambda}  k_{2\sigma} k_{4\lambda}
\Big[ - k_{1} \cdot k_{2} G_B(u_1, u_2)
 \nonumber\\
 &&- k_{1} \cdot k_{3} G_B(u_1, u_3) - k_{1}\cdot k_{4} G_B(u_1, u_4) - k_{2} \cdot k_{3} G_B(u_2, u_3) - k_{2} \cdot k_{4} G_B(u_2, u_4) - k_{3} \cdot k_{4} G_B(u_3, u_4) \Big]^{-2}\,.
 \label{Gammau}
\end{eqnarray}

Finally, we need to integrate over the proper time variables $u_k$ which define the position of the interaction points on the worldline. Using the rotational invariance of the worldline loop~\cite{Schubert:2001he}, it is convenient to fix $u_1 = 0$ (which is equivalent then to $u_1=1$). There are six possible orderings~\cite{Kodaira:1979pa} of the proper time variables $u_k$, which can be split into the two classes  shown in Fig. \ref{fig:box_ordering}a and \ref{fig:box_ordering}b.
\begin{figure}[htb]
 \begin{center}
 \includegraphics[width=140mm]{box_ordering}
 \end{center}
 \caption{\label{fig:box_ordering}Two distinct topologies corresponding to the ordering of the proper time coordinates $u_k$.}
 \end{figure}
 
The first class of diagrams in Fig. \ref{fig:box_ordering}a corresponds to configurations  where there are no gluon insertions between the two electromagnetic currents. There are four such possible orderings: $u_3 > u_2 > u_4$, $u_3 > u_4 > u_2$, $u_2 > u_4 > u_3$, $u_4 > u_2 > u_3$.
By explicit calculation of the integrals over $u_k$ presented below, one finds that all these orderings yield the same contribution with an infrared anomaly pole $\sim \frac{1}{(k_2 + k_4)^2}$. Indeed, as discussed at length in \cite{Usyukina:1992jd}, such graphs have the  generic structure $\frac{1}{(k_2+k_4)^2}\,\frac{1}{(k_1+k_2)^2}\rightarrow \frac{1}{t}\,\frac{1}{Q^2+2xP\cdot q}$, where recall $k_2+k_4=l$,  
$k_1^2= Q^2$ and we define $k_2 = 2\,x P$.

The second class of diagrams in Fig. \ref{fig:box_ordering}b corresponds to a gluon insertion between the two electromagnetic currents (corresponding to a ``cat's eye" topology) and have two possible orderings: $u_2 > u_3 > u_4$ and $u_4 > u_3 > u_2$. These diagrams do not have an infrared pole ! This is because in this case the diagrams have the generic structure $\frac{1}{(k_1+k_4)^2}\,\frac{1}{(k_2+k_3)^2}\rightarrow \frac{1}{Q^2+2xP\cdot q}\,\frac{1}{Q^2+2xP\cdot q}$  and therefore do not have an infrared pole in the forward limit $l_\mu\to 0$. Note that this class of diagrams is suppressed by a factor of $1/Q^2$ relative to those given by Fig. \ref{fig:box_ordering}a and therefore does not contribute  in the Bjorken limit. 


Since all four orderings in Fig. \ref{fig:box_ordering}a give the same result, we will fix the ordering of the proper times variables  (multiplying the expression by a factor of $4$) in Eq. (\ref{Gammau}) as
\begin{eqnarray}
&&\Gamma^{\mu\nu\alpha\beta}_A[k_1, k_3, k_2, k_4]\Big|_{Q^2\to \infty}
= 4 \frac{g^2e^2 e_f^2 }{4\pi^2}\epsilon^{\mu\nu\eta}_{\ \ \ \ \kappa} (k_{1\eta} - k_{3\eta})(2\pi)^4\delta^{(4)}(\sum^4_{i=1} k_i)
 \int^1_0 du_3 \int^{u_3}_0 du_4 \int^{u_4}_0 du_2 
 \nonumber\\
&&\times \Big[ - G_B(0, u_2) + G_B(u_2, u_4) - G_B(u_4, 0)
 \Big] ( k^\kappa_{2} + k^\kappa_{4} ) \epsilon^{\alpha\beta\sigma\lambda}  k_{2\sigma} k_{4\lambda}
\Big[ - k_{1} \cdot k_{2} G_B(0, u_2)
 \nonumber\\
 &&- k_{1} \cdot k_{3} G_B(0, u_3) - k_{1}\cdot k_{4} G_B(0, u_4) - k_{2} \cdot k_{3} G_B(u_2, u_3) - k_{2} \cdot k_{4} G_B(u_2, u_4) - k_{3} \cdot k_{4} G_B(u_3, u_4) \Big]^{-2}\,.
 \label{gammaordered}
\end{eqnarray}

We now introduce the change of variables $u_i\rightarrow a_i$, where the latter can be identified as the standard Feynman parameters~\cite{Strassler:1992zr},
\begin{eqnarray}
a_1 = 1 - u_3;\ \ \ a_2 = u_3 - u_4;\ \ \ a_3 = u_4 - u_2;\ \ \ a_4 = u_2\,,
\end{eqnarray}
and rewrite Eq. (\ref{gammaordered}) as
\begin{eqnarray}
&&\Gamma^{\mu\nu\alpha\beta}_A[k_1, k_3, k_2, k_4]\Big|_{Q^2\to \infty}
= 4 \frac{g^2e^2 e_f^2 }{4\pi^2}\epsilon^{\mu\nu\eta}_{\ \ \ \ \kappa} (k_{1\eta} - k_{3\eta})(2\pi)^4\delta^{(4)}(\sum^4_{i=1} k_i)
 \prod^4_{k=1}\int^1_0 da_k ~\delta(1 - \sum^4_{j=1}a_j)
 \nonumber\\
&&\times \Big[ - 2a_2a_4 \Big] ( k^\kappa_{2} + k^\kappa_{4} ) \epsilon^{\alpha\beta\sigma\lambda}  k_{2\sigma} k_{4\lambda}
\Big[ (k_{2} + k_4)^2  a_2a_4 + ( k_{1} + k_2)^2 a_1 a_3 + k^2_{1}  a_1a_4 + k^2_{3} a_1 a_2 \Big]^{-2}\,.
\end{eqnarray}
We have neglected terms $\sim a_1$ in the numerator of this equation; as mentioned previously, these terms, in the Bjorken limit, are  of $O(1/Q^4)$ and are therefore suppressed. To see this more clearly, following \cite{Usyukina:1992jd}, we introduce a further change of variables,
\begin{eqnarray}
a_1 = \alpha;\ \ \ a_i = (1-\alpha)\beta_i,\ \ \ i=2,3,4\,,
\end{eqnarray}
and rewrite the equation as
\begin{eqnarray}
&&\Gamma^{\mu\nu\alpha\beta}_A[k_1, k_3, k_2, k_4]\Big|_{Q^2\to \infty}
= 4 \frac{g^2e^2 e_f^2 }{4\pi^2}\epsilon^{\mu\nu\eta}_{\ \ \ \ \kappa} (k_{1\eta} - k_{3\eta})(2\pi)^4\delta^{(4)}(\sum^4_{i=1} k_i)
 \int^1_0 d\alpha ~\prod^4_{k=2}\int^1_0 d\beta_k ~\delta(1 - \sum_{k=2,3,4}\beta_k ) 
 \nonumber\\
&&\times  \Big[ - 2 \beta_2\beta_4 \Big] ( k^\kappa_{2} + k^\kappa_{4} ) \epsilon^{\alpha\beta\sigma\lambda}  k_{2\sigma} k_{4\lambda}
\Big[ (k_{2} + k_4)^2  (1 - \alpha)\beta_2\beta_4 + ( k_{1} + k_2)^2 \alpha \beta_3 + k^2_{1}  \alpha \beta_4 + k^2_{3} \alpha \beta_2 \Big]^{-2}\,,
\end{eqnarray}

The integral over $\alpha$ can be performed easily, giving
\begin{eqnarray}
\Gamma^{\mu\nu\alpha\beta}_A[k_1, k_3, k_2, k_4]\Big|_{Q^2\to \infty}
&=& - 2\frac{g^2e^2 e_f^2 }{\pi^2}\epsilon^{\mu\nu\eta}_{\ \ \ \ \kappa} (k_{1\eta} - k_{3\eta}) (2\pi)^4\delta^{(4)}(\sum^4_{i=1} k_i)
 ~\frac{( k^\kappa_{2} + k^\kappa_{4} ) \epsilon^{\alpha\beta\sigma\lambda}  k_{2\sigma} k_{4\lambda} }{(k_{2} + k_4)^2  }  
 \nonumber\\
&\times&  \prod^4_{k=2}\int^1_0 d\beta_k ~\delta(1 - \sum_{k=2,3,4}\beta_k ) \frac{1}{ ( k_{1} + k_2)^2 \beta_3 + k^2_{1} \beta_4 + k^2_{3} \beta_2}\,.
\label{Gammawithinfpole}
\end{eqnarray}
Note that the first line on the r.h.s of this equation contains our result (computed explicitly in Appendix C) for the triangle diagram in Eq.~(\ref{eq:VVA1}). 
Specifically, we see that Eq. (\ref{Gammawithinfpole}) has the infrared anomaly pole $\frac{1}{(k_{2} + k_4)^2}\equiv\frac{1}{t}$, which doesn't depend on the integration over $\beta_k$. Hence it is safe to take the forward limit in the integrals over $\beta_k$, and in particular, the forward limit relation between the  momenta of the incoming virtual photons: $k^2_3 = k^2_1$. One can then simplify the equation to read,
\begin{eqnarray}
\Gamma^{\mu\nu\alpha\beta}_A[k_1, k_3, k_2, k_4]\Big|_{Q^2\to \infty}
&=& - 2\frac{g^2e^2 e_f^2 }{\pi^2}\epsilon^{\mu\nu\eta}_{\ \ \ \ \kappa} (k_{1\eta} - k_{3\eta}) (2\pi)^4\delta^{(4)}(\sum^4_{i=1} k_i)
 ~\frac{( k^\kappa_{2} + k^\kappa_{4} ) \epsilon^{\alpha\beta\sigma\lambda}  k_{2\sigma} k_{4\lambda} }{(k_{2} + k_4)^2  }  
 \nonumber\\
&\times&  \prod^4_{k=2}\int^1_0 d\beta_k ~\delta(1 - \sum_{k=2,3,4}\beta_k ) \frac{1}{ 2 k_1 \cdot k_2 \beta_3   + k^2_{1} }\,.
\end{eqnarray}
Next, one can straightforwardly perform the integration over the $\beta_k$ variables, which gives
\begin{eqnarray}
\Gamma^{\mu\nu\alpha\beta}_A[k_1, k_3, k_2, k_4]\Big|_{Q^2\to \infty}
&=& - 2\frac{g^2e^2 e_f^2 }{\pi^2}\epsilon^{\mu\nu\eta}_{\ \ \ \ \kappa} (k_{1\eta} - k_{3\eta}) (2\pi)^4\delta^{(4)}(\sum^4_{i=1} k_i)
 ~\frac{( k^\kappa_{2} + k^\kappa_{4} ) \epsilon^{\alpha\beta\sigma\lambda}  k_{2\sigma} k_{4\lambda} }{(k_{2} + k_4)^2  }  
 \nonumber\\
&\times&  \frac{1}{2 k_1 \cdot k_2}\Big[ \big( 1  + \frac{k^2_{1}}{2 k_1 \cdot k_2} \big)\ln\Big[ \frac{2 k_1 \cdot k_2 + k^2_{1}}{k^2_{1}}\Big]  - 1 \Big]\,,
\label{Gammaaftb}
\end{eqnarray}
where we used the identity
\begin{eqnarray}
\int^1_0 \frac{dx}{ax + b} = \frac{1}{a}\ln\Big[\frac{a+b}{b}\Big]\,,
\end{eqnarray}
employed in standard computations of the box diagram~\cite{tHooft:1978jhc}.

Now substituting the above expression into the box diagram (Eq. (\ref{Hadgamma})), we can write the antisymmetric piece of the polarization tensor in the Bjorken limit as 
\begin{eqnarray}
\label{Gafin}
&&\Gamma^{\mu\nu}_A[k_1, k_3]\Big|_{Q^2\to \infty} = -\frac{g^2e^2 e_f^2 }{\pi^2}\epsilon^{\mu\nu\eta}_{\ \ \ \ \kappa} (k_{1\eta} - k_{3\eta}) ~ \frac{ k^\kappa_{2} + k^\kappa_{4} }{(k_{2} + k_4)^2  }
\nonumber \\
&&\times \int \frac{d^4k_2}{(2\pi)^4} \int \frac{d^4k_4}{(2\pi)^4}~
(2\pi)^4 \delta^{(4)}(\sum^4_{i=1} k_i)
~  \frac{1}{2 k_1 \cdot k_2}\Big[ \big( 1  + \frac{k^2_{1}}{2 k_1 \cdot k_2} \big)\ln\Big[ \frac{2 k_1 \cdot k_2 + k^2_{1}}{k^2_{1}}\Big]  - 1 \Big]~ {\rm Tr_c}(F_{\alpha\beta}(k_2) \tilde{F}^{\alpha\beta}(k_4))\,.
\end{eqnarray}
In arriving at this expression, we have generalized the expression in terms of derivatives of the background fields to express it in terms of the field strength tensors $F^{\mu\nu}$ and its dual ${\tilde F}^{\mu\nu}$, where recall $\tilde{F}_{\mu\nu} = \frac{1}{2}\,\epsilon_{\mu\nu\rho\sigma}F^{\rho\sigma}$.

Finally, substituting Eq. (\ref{Gafin}) into the antisymmetric piece of the hadron tensor (Eq.~(\ref{WModFu})), we obtain
\begin{eqnarray}
i\tilde{W}_f^{\mu\nu}(q, P, S) \Big|_{Q^2\to \infty}
&=& - \frac{g^2 e_f^2 }{\pi^3}\,{\rm Im}~ \epsilon^{\mu\nu\eta}_{\ \ \ \ \kappa} q_{\eta} ~ 
~ \int \frac{d^4k}{(2\pi)^4}~
~  \frac{1}{ 2 q \cdot k}\Big[ \Big( 1  + \frac{q^2}{2 q \cdot k} \Big)\ln\Big[ \frac{2 q \cdot k + q^2}{q^2}\Big]  - 1 \Big]
\nonumber\\
&\times& \int d^4 z\, e^{-ik z } \lim_{l_\kappa \to 0}\frac{ l^\kappa }{l^2  }\langle P^\prime,S| {\rm Tr_c} F_{\alpha\beta}(z) \tilde{F}^{\alpha\beta}(0) |P,S\rangle\,,
\label{WantisimA}
\end{eqnarray}
where ${\tilde W}^{\mu\nu}(q, P, S)=\sum_f\tilde{W}_f^{\mu\nu}(q, P, S)$ and  ${P^\prime}^\kappa-P^\kappa = l^\kappa$. 
In writing this result, we used Eq. (\ref{fourbackground}), integrated over intermediate momenta and coordinates, and performed an analytical continuation of the expression to Minkowski space-time.

Following  \cite{Leader:2001gr}, we will make a high energy approximation and write $k^\mu \approx (k\cdot n) P^\mu$, where $n$ is a dimensionful vector such that $n^2 = 0$ and $n\cdot P = 1$. As a result, we can rewrite the equation as
\begin{eqnarray}
i\tilde{W}_f^{\mu\nu}(q, P, S) \Big|_{Q^2\to \infty}
&=& - \frac{g^2 e_f^2 }{\pi^3}\,{\rm Im} \frac{1}{ 2 P \cdot q} \epsilon^{\mu\nu\eta}_{\ \ \ \ \kappa} q_{\eta} ~
 \int^1_0 dx 
~  \frac{1}{x}\Big[ \Big( 1  - \frac{x_B}{x} \Big)\ln\Big[ \frac{x - x_B}{ - x_B} \Big]  - 1 \Big]
\nonumber\\
&\times& \int \frac{d^4k}{(2\pi)^4}~\delta(x - k\cdot n) \lim_{l_\kappa \to 0}\frac{ l^\kappa }{l^2  } \int d^4 z\, e^{-ik z } \langle P^\prime,S| {\rm Tr_c} F_{\alpha\beta}(z) \tilde{F}^{\alpha\beta}(0) |P,S\rangle\,. 
\end{eqnarray}
 The imaginary part of the expression on the r.h.s is obtained from the identity $\ln(x) = \ln(|x|) + i\pi$, which also requires that we impose $x\geq x_B$. Hence,
\begin{eqnarray}
i\tilde{W}_f^{\mu\nu}(q, P, S) \Big|_{Q^2\to \infty}
&=& \frac{g^2 e_f^2 }{\pi^2} \frac{1}{ 2 P \cdot q} \epsilon^{\mu\nu\eta}_{\ \ \ \ \kappa} q_{\eta} ~ 
  \int^1_{x_B} \frac{dx}{x} 
~ \Big( 1  - \frac{x_B}{x} \Big)
\nonumber\\
&\times& \int \frac{d^4k}{(2\pi)^4}~\delta(x - k\cdot n) \lim_{l_\kappa \to 0}\frac{ l^\kappa }{l^2  }\,\int d^4 z e^{-ik z } \langle P^\prime,S| {\rm Tr_c} F_{\alpha\beta}(z) \tilde{F}^{\alpha\beta}(0) |P,S\rangle\,.
\end{eqnarray}
Integrating over $k^\mu$ further yields,
\begin{eqnarray}
i\tilde{W}_f^{\mu\nu}(q, P, S) \Big|_{Q^2\to \infty}
= \frac{4}{ P \cdot q} \epsilon^{\mu\nu\eta}_{\ \ \ \ \kappa} q_{\eta} ~
 \int^1_{x_B} \frac{dx}{x} 
~ \Big( 1  - \frac{x_B}{x} \Big) \int \frac{d\xi}{2\pi} e^{-i\xi x } \lim_{l_\kappa \to 0}\frac{ l^\kappa }{l^2  } \langle P^\prime,S| \frac{\alpha_s e_f^2}{2\pi} {\rm Tr_c} F_{\alpha\beta}(\xi n) 
\tilde{F}^{\alpha\beta}(0) |P,S\rangle\,,
\nonumber\\
\label{WAresBj}
\end{eqnarray}
where the r.h.s is rewritten in terms of the matrix element of the nonlocal operator 
\begin{eqnarray}
\int \frac{d\xi}{2\pi} e^{-i\xi x }\lim_{l_\kappa \to 0}\frac{ l^\kappa }{l^2  }  \langle P^\prime,S| \frac{\alpha_s e_f^2}{2\pi} {\rm Tr_c} F_{\alpha\beta}(\xi n) \tilde{F}^{\alpha\beta}(0) |P,S\rangle\,.
\end{eqnarray}
Comparing Eq. (\ref{WAresBj}) with the general tensorial decomposition of the antisymmetric part of the hadron tensor in Eq.~(\ref{WA}), allows us to extract our final result for $g_1(x_B, Q^2)$ in the Bjorken limit:
\begin{eqnarray}
\label{eq:g1-Bj-final}
 S^\mu g_1(x_B, Q^2)\Big|_{Q^2\to \infty} = \sum_f e_f^2 \frac{\alpha_s}{i\pi M_N}
 \int^1_{x_B} \frac{dx}{x} 
~ \Big( 1  - \frac{x_B}{x} \Big) \int \frac{d\xi}{2\pi} e^{-i\xi x }\lim_{l_\mu \to 0} \frac{ l^\mu }{l^2  }\langle P^\prime,S| {\rm Tr_c} F_{\alpha\beta}(\xi n) \tilde{F}^{\alpha\beta}(0) |P,S\rangle\,.
\label{resultBj}
\end{eqnarray}

Finally, as noted in the introduction, the first moment of $g_1(x_B,Q^2)$ is simply related to $\Sigma(Q^2)$ and is given by 
\begin{eqnarray}
 S^\mu\int^1_0 dx_B \,g_1(x_B, Q^2)\Big|_{Q^2\to \infty} = \sum_f e_f^2 \frac{\alpha_S}{2 i\pi M_N}\,
 \lim_{l_\mu \to 0} \frac{ l^\mu }{l^2  } \langle P^\prime ,S| {\rm Tr_c} F_{\alpha\beta}(0) \tilde{F}^{\alpha\beta}(0) |P,S\rangle\,.
\end{eqnarray}
Note that this expression is the first term on the r.h.s of Eq.~(5.19) in Ref.~\cite{Jaffe:1989jz}.  Clearly, this differs from Eq.~(\ref{eq:sigma-anomaly}) because, as discussed in the introduction, the contribution from the pseudoscalar sector is not included here.

\subsection{Discussion of Bjorken limit result in Eq.~(\ref{eq:g1-Bj-final})}
\label{subsec:Bj-discussion}

The result for $g_1(x_B, Q^2)$ is one of the principal results of this paper. It shows that the box diagram for polarized DIS in the Bjorken limit is dominated by the triangle anomaly and if not regulated appropriately will diverge in the forward limit $l^\kappa\rightarrow 0$. Though nonlocal, the nonlocality must be interpreted as a smearing of the operator corresponding to the bare topological charge density $Q_B = \frac{\alpha_S}{8\pi} {\rm Tr}\left(F {\tilde F}\right)$. This is because our derivation parallels exactly the derivation (worked out in Appendix~\ref{sec:triangle}) of the triangle anomaly in the worldline formalism. Indeed, as noted earlier, the first line of Eq.~(\ref{Gammawithinfpole}) already contains the infrared pole of the anomaly, with the terms on the second line only giving finite contributions. 

As we will discuss further in Papers II\,\& III, this operator $Q_B$ undergoes renormalization with evolution  and mixes, via the anomaly, with the $t$-channel exchange of a massless isosinglet pseudoscalar $\eta_0$ to generate the massive $\eta^\prime$ meson. Thus as emphasized by Veneziano~\cite{Veneziano:1989ei}, how one recovers a finite result for $g_1(x_B,Q^2)$ in the forward limit $l^\kappa\rightarrow 0$ is deeply tied to the resolution of the $U_A(1)$ problem in QCD. This is independent of whether the underlying dynamical mechanism is  due to instantons or other nonperturbative phenomena and a consequence of anomalous chiral Ward identities~\cite{Shore:2007yn}.

With regard to the underlying dynamical mechanism, we should emphasize that the computation of the matrix element on the r.h.s of the above expression is nontrivial. Firstly, if we  naively take the forward limit, the matrix element vanishes. One way to properly interpret this matrix element (along with the factor $l^\mu/l^2$) is to consider its path integral realization as a convolution of the operator and the density matrix of states. As is well known, the  anomaly arises from the measure of the path integral~\cite{Fujikawa:1979ay}. In the worldline formalism, this contribution of the measure can be reexpressed as the imaginary part of the one loop effective action (or equivalently, the phase of the fermion determinant). Thus the  effect of the anomaly is distinct from that of the real part of the effective action, which is, for instance, responsible for the equations of motion ~\cite{Mueller:2017arw,Mueller:2019gjj}. 

In addition to the anomaly contribution,  the imaginary part of the effective action contains a Wess-Zumino-Witten term for the $\eta_0$, as noted in the introduction, that will cancel the pole of the anomaly in the computation of the path integral~\cite{DHoker:1995aat}. This key point will be discussed further in Paper II. Our perspective is  similar to that of Shore, Narison and  Veneziano~\cite{Narison:1994hv,Shore:2007yn, Narison:1998aq}, who employed the Wess-Zumino effective action~\cite{Wess:1971yu} in their studies to obtain the same results for the first moment of $g_1(x_B,Q^2)$.

To discuss further the implications of our result, observe  that the matrix element in Eq.~(\ref{eq:g1-Bj-final}) has a distinct tensorial structure from that in the polarized gluon distribution\footnote{We omit writing the gauge link between $F^{+\mu}$ and ${\tilde F}^+_{\ \mu}$ required to ensure gauge invariance.}, which is defined as~\cite{Kodaira:1998jn}
\begin{eqnarray}
\label{eq:DeltaGx_B}
\Delta G(x_B,Q^2) &=& \frac{2i}{x_B}\int \frac{d\xi}{2\pi}\,e^{-i\xi x_B} \langle P,S| {\rm Tr_c} n_\alpha F^{\alpha\mu}(\xi n) n^\beta \tilde{F}_{\beta\mu}(0) |P,S\rangle\nonumber \\
&\equiv& \frac{2i}{x_B (P^+)^2}\int \frac{d\xi}{2\pi}\,e^{-i\xi x_B} \langle P,S| {\rm Tr_c} F^{+\mu}(\xi n) \tilde{F}^+_{\ \mu}(0) |P,S\rangle\,.
\end{eqnarray}
More importantly, as we have argued, $g_1(x_B,Q^2)$ in Eq.~(\ref{eq:g1-Bj-final}) is dominated entirely by the anomaly contribution from the imaginary part of the worldline effective action in the path integral for the matrix element. It is therefore unclear how relate it to the r.h.s of the above expression\footnote{We thank Bob Jaffe for bringing his Varenna lectures (unpublished) to our attention, where this point is discussed at some length.} as is often done in perturbative computations in the literature~\cite{Zijlstra:1993sh,Ball:1995ye,Blumlein:1996hb}. Further clarity\footnote{We thank Y. Hatta for an interesting discussion on this point.} on this important issue can be obtained from the explicit computation of the matrix elements for $g_1(x_B,Q^2)$ and $\Delta G(x_B,Q^2)$.

Our results are consistent with perturbative computations of the renormalization group evolution of $\Sigma(Q^2)$ with $Q^2$. 
In pQCD computations, the RG evolution of $\Sigma(Q^2)$ can mix that of the first moment $\Delta G(Q^2)$ of Eq.~(\ref{eq:DeltaGx_B}), the matrix element of the other isosinglet twist two operator in polarized DIS. Their combined evolution is described by a two-by-two matrix of splitting functions of polarized quarks splitting into softer polarized quarks or gluons and likewise, polarized gluons splitting into softer polarized gluons or quarks. These splitting functions were  computed to leading order (LO) in \cite{Ahmed:1976ee, Altarelli:1977zs,Kodaira:1979pa}, to next-to-leading order (NLO) in \cite{Mertig:1995ny, Vogelsang:1995vh,Vogelsang:1996im} and to next-to-next-to-leading order (NNLO) in 
\cite{Moch:2014sna,Moch:2015usa}. As shown already for the leading order splitting functions in \cite{Altarelli:1990jp}, but remarkably also for the values of the splitting functions computed to NNLO accuracy~\cite{deFlorian:2019egz},  $\Sigma(Q^2)$ does not mix with $\Delta G(Q^2)$. This is precisely what we would expect since $\Sigma(Q^2)$ has a distinct topological structure determined by the chiral anomaly alone. On the other hand, while $\Delta G$ does not influence the evolution of $\Sigma(Q^2)$, the converse needn't be true, as also observed in the pQCD computations. We would argue that this is because $\Delta G$ is sensitive to fermion  zero modes that are influenced by the anomaly. 

\section{The triangle anomaly in the Regge limit of the box diagram}
\label{sec:Regge}
We  now turn our attention to the computation of the box diagram in the other interesting asymptotics of QCD, the Regge limit of $x_B\ll 1$ for fixed $Q^2\gg \Lambda_{\rm QCD}^2$.  In this regime, the quarks in the box diagram suffer a nearly instantaneous shock wave interaction with  the background gluons, which corresponds to taking the limit $u_2\rightarrow u_4$ in Eq.~(\ref{GammawithIinu}), given by Eq.~(\ref{GammawithIRj}). As we shall see, the derivation closely parallels that of our discussion of the Bjorken limit, with some small but important differences that we will draw the reader's attention to. 

\subsection{Worldline computation of the box diagram in the Regge limit}
\label{subsec:Rj-worldline}

A cursory examination reveals that taking the limit $u_2\rightarrow u_4$ in the box diagram is very similar to taking  $u_1\rightarrow u_3$ in the Bjorken limit. Thus one can obtain the sum 
of the coefficients in Eq. (\ref{GammawithIRj}) by borrowing Eq. (\ref{suminBj}) from the previous section and making the following substitutions:  $\mu\leftrightarrow\alpha$, $\nu\leftrightarrow\beta$, $k_1\leftrightarrow k_2$, $k_3\leftrightarrow k_4$, $u_1 \leftrightarrow u_2$, $u_3\leftrightarrow u_4$:
\begin{eqnarray}
\label{eq:Regge-tensor-sum}
&&\Big(\sum^9_{n=1}\mathcal{I}^{\mu\nu\alpha\beta}_{n;(u_1,u_2,u_3,u_4)}[k_1, k_3, k_2, k_4]\Big)\Big|_{u_2 = u_4} = \epsilon^{\alpha\beta\eta}_{\ \ \ \ \kappa} (k_{2\eta} - k_{4\eta})
\nonumber\\
&&\times \Big( k_1 \cdot k_3 \Big[ \dot{G}^2_B(u_2, u_3) - \dot{G}_B(u_2, u_3) \dot{G}_B(u_1, u_3) + \dot{G}_B(u_2, u_1) \big(\dot{G}_B(u_2, u_1) + \dot{G}_B(u_1, u_3) + \dot{G}_B(u_3, u_2) \big)  \Big] \epsilon^{\kappa\mu\nu\sigma} k_{1\sigma}
\nonumber\\
&&+ k_1\cdot k_3 \Big[ -\dot{G}^2_B(u_2, u_1) - \dot{G}_B(u_2, u_1) \dot{G}_B(u_1, u_3) + \dot{G}_B(u_2, u_3) \big(\dot{G}_B(u_2, u_1) + \dot{G}_B(u_1, u_3) + \dot{G}_B(u_3, u_2)\big) \Big] \epsilon^{\kappa\mu\nu\sigma} k_{3\sigma}
\nonumber\\
&&+ \Big[ - \dot{G}_B(u_2, u_3) \dot{G}_B(u_1, u_3) +  \dot{G}^2_B(u_2, u_3)
+ \dot{G}_B(u_2, u_1) \big(\dot{G}_B(u_2, u_1) + \dot{G}_B(u_1, u_3) + \dot{G}_B(u_3, u_2)\big) \Big] \epsilon^{\kappa\mu\sigma\lambda} k^\nu_{1} k_{1\sigma} k_{3\lambda}
\nonumber\\
&&+ \dot{G}^2_B(u_2, u_3) \epsilon^{\kappa\mu\sigma\lambda} k^\nu_{3} k_{1\sigma} k_{3\lambda}
- \dot{G}^2_B(u_2, u_1) \epsilon^{\kappa\nu\sigma\lambda} k^\mu_{1} k_{1\sigma} k_{3\lambda}
\nonumber\\
&&+ \Big[ - \dot{G}_B(u_2, u_1) \dot{G}_B(u_1, u_3) - \dot{G}^2_B(u_2, u_1)
+ \dot{G}_B(u_2, u_3) \big(\dot{G}_B(u_2, u_1) + \dot{G}_B(u_1, u_3) + \dot{G}_B(u_3, u_2)\big) \Big] \epsilon^{\kappa\nu\sigma\lambda} k^\mu_{3} k_{1\sigma} k_{3\lambda}
\nonumber\\
&&- \epsilon^{\mu\lambda\nu\sigma}  k^\kappa_{1} k_{1\lambda} k_{3\sigma} - \epsilon^{\mu\lambda\nu\sigma} k^\kappa_{3} k_{1\lambda} k_{3\sigma}\Big)\,.
\end{eqnarray}
In analogy to the previous derivation, this expression can be simplified by  employing the cyclic identity for the Levi-Civita tensor (in footnote~\ref{cyclic}), which gives,
\begin{eqnarray}
\label{eq:LC_R1}
- \epsilon^{\mu\lambda\nu\sigma}  k^\kappa_{1} k_{1\lambda} k_{3\sigma} = - k_{1}\cdot k_{3} \epsilon^{\kappa\mu\nu\sigma} k_{1\sigma} - \epsilon^{\kappa\mu\sigma\lambda} k^\nu_{1} k_{1\sigma} k_{3\lambda} + \epsilon^{\kappa\nu\sigma\lambda} k^\mu_{1} k_{1\sigma} k_{3\lambda} + k^2_{1} \epsilon^{\nu\sigma\kappa\mu} k_{3\sigma}\,,
\end{eqnarray}
and
\begin{eqnarray}
\label{eq:LC_R2}
&&- \epsilon^{\mu\lambda\nu\sigma} k^\kappa_{3} k_{1\lambda} k_{3\sigma} = k_{1}\cdot k_{3} \epsilon^{\kappa\mu\nu\sigma} k_{3\sigma} + \epsilon^{\kappa\nu\sigma\lambda} k^\mu_{3} k_{1\sigma} k_{3\lambda} - \epsilon^{\kappa\mu\sigma\lambda} k^\nu_{3} k_{1\sigma} k_{3\lambda} - k^2_{3} \epsilon^{\kappa\mu\nu\sigma} k_{1\sigma}\,.
\end{eqnarray}
Note that unlike  Eq.~(\ref{Levi-Civita}), where the term corresponding to the virtuality of the background gluons was set to be $k_2^2=k_4^2=0$ in the Bjorken limit, here the analogous terms give $k_1^3= k_3^2=Q^2\neq 0$.
With Eqs.~(\ref{eq:LC_R1}) and (\ref{eq:LC_R2}), we can rewrite Eq.~(\ref{eq:Regge-tensor-sum}) as
\begin{eqnarray}
&&\Big(\sum^9_{n=1}\mathcal{I}^{\mu\nu\alpha\beta}_{n;(u_1,u_2,u_3,u_4)}[k_1, k_3, k_2, k_4]\Big)\Big|_{u_2 = u_4} = \epsilon^{\alpha\beta\eta}_{\ \ \ \ \kappa} (k_{2\eta} - k_{4\eta})
\nonumber\\
&&\times \Big( \Big[ 1 - \dot{G}^2_B(u_2, u_3) + \dot{G}_B(u_2, u_3) \dot{G}_B(u_1, u_3) - \dot{G}_B(u_2, u_1) \big(\dot{G}_B(u_2, u_1) + \dot{G}_B(u_1, u_3) + \dot{G}_B(u_3, u_2) \big)  \Big] k^\kappa_{1} \epsilon^{\mu\sigma\lambda\nu} k_{1\sigma} k_{3\lambda}
\nonumber\\
&&+ \Big[ -1 + \dot{G}^2_B(u_2, u_1) + \dot{G}_B(u_2, u_1) \dot{G}_B(u_1, u_3) - \dot{G}_B(u_2, u_3) \big(\dot{G}_B(u_2, u_1) + \dot{G}_B(u_1, u_3) + \dot{G}_B(u_3, u_2)\big) \Big]k^\kappa_{3} \epsilon^{\nu\sigma\lambda\mu} k_{1\sigma} k_{3\lambda}
\nonumber\\
&&+ \Big[ - \dot{G}^2_B(u_2, u_1) - \dot{G}_B(u_2, u_1) \dot{G}_B(u_1, u_3) +  \dot{G}_B(u_2, u_3) \dot{G}_B(u_2, u_1) + \dot{G}_B(u_2, u_3) \dot{G}_B(u_1, u_3) \Big] \epsilon^{\kappa\mu\sigma\lambda} k^\nu_{3} k_{1\sigma} k_{3\lambda}
 \nonumber\\
 &&+ \Big[ \dot{G}^2_B(u_2, u_3) - \dot{G}_B(u_2, u_3) \dot{G}_B(u_1, u_3) 
 + \dot{G}_B(u_2, u_1) \dot{G}_B(u_1, u_3) + \dot{G}_B(u_2, u_1) \dot{G}_B(u_3, u_2) \Big] \epsilon^{\kappa\nu\sigma\lambda} k^\mu_{1} k_{1\sigma} k_{3\lambda} 
 \nonumber\\
 &&+\Big[ \dot{G}^2_B(u_2, u_3) - \dot{G}_B(u_2, u_3) \dot{G}_B(u_1, u_3) + \dot{G}_B(u_2, u_1) \big(\dot{G}_B(u_2, u_1) + \dot{G}_B(u_1, u_3) + \dot{G}_B(u_3, u_2) \big)  \Big] k^2_{1} \epsilon^{\nu\sigma\kappa\mu} k_{3\sigma}
 \nonumber\\
&&+ \Big[ \dot{G}^2_B(u_2, u_1) + \dot{G}_B(u_2, u_1) \dot{G}_B(u_1, u_3) - \dot{G}_B(u_2, u_3) \big(\dot{G}_B(u_2, u_1) + \dot{G}_B(u_1, u_3) + \dot{G}_B(u_3, u_2)\big) \Big] k^2_{3} \epsilon^{\mu\kappa\nu\sigma} k_{1\sigma}
\end{eqnarray}


From gauge invariance, the terms proportional to $\sim k^\mu_1$ and $\sim k^\nu_3$ do not  contribute to the final answer. Further, the terms proportional to  $\sim k^2_1$ and $\sim k^2_3$ will be proportional to $ \frac{Q^2}{2P\cdot q} \equiv x_B$ in the forward limit and are  suppressed, relative to the other terms, in Regge asymptotics.  The leading contributions are therefore given by the first two lines in the previous equation, which can be written in a greatly simplified form as 
\begin{eqnarray}
&&\Big(\sum^9_{n=1}\mathcal{I}^{\mu\nu\alpha\beta}_{n;(u_1,u_2,u_3,u_4)}[k_1, k_3, k_2, k_4]\Big)\Big|_{u_2 = u_4} = - \frac{1}{2} \epsilon^{\alpha\beta\eta}_{\ \ \ \ \kappa} (k_{2\eta} - k_{4\eta})
\nonumber \\
&&\times \Big[ -2 + \dot{G}^2_B(u_2, u_1) - \dot{G}^2_B(u_1, u_3) + \dot{G}^2_B(u_3, u_2)
+ \big(\dot{G}_B(u_2, u_1) + \dot{G}_B(u_1, u_3) + \dot{G}_B(u_3, u_2) \big)^2 \Big] 
(k^\kappa_{1} + k^\kappa_{3}) \epsilon^{\mu\nu\sigma\lambda} k_{1\sigma} k_{3\lambda}\,.
\nonumber\\
\end{eqnarray}

Using Eqs. (\ref{gfprod}) and (\ref{b2dot}) to further simplify the expression in the brackets, we obtain 
\begin{eqnarray}
\Big(\sum^9_{n=1}\mathcal{I}^{\mu\nu\alpha\beta}_{n;(u_1,u_2,u_3,u_4)}[k_1, k_3, k_2, k_4]\Big)\Big|_{u_2 = u_4} 
&=& - 2 \epsilon^{\alpha\beta\eta}_{\ \ \ \ \kappa} (k_{2\eta} - k_{4\eta})
 \Big[ - G_B(u_2, u_1) + G_B(u_1, u_3) - G_B(u_3, u_2) \Big] \nonumber \\
 &\times& (k^\kappa_{1} + k^\kappa_{3}) \epsilon^{\mu\nu\sigma\lambda} k_{1\sigma} k_{3\lambda}\,.
\end{eqnarray}
Finally,  we can reexpress this result as\footnote{We employ here energy-momentum conservation ($k_1^\kappa+k_3^\kappa = - k_2^\kappa - k_4^\kappa$), and the identity $\epsilon^{\mu\nu\sigma\lambda} k_{1\sigma}k_{3\sigma} = \frac{1}{2}\epsilon^{\mu\nu\sigma\lambda} (k_1-k_3)_\sigma (k_1+k_3)_\lambda$  (likewise, $\epsilon^{\alpha\beta\sigma\lambda} (k_2-k_4)_\sigma (k_2+k_4)_\lambda = 2 \epsilon^{\alpha\beta\sigma\lambda} k_{2\sigma} k_{4\lambda}$).} 
\begin{eqnarray}
\Big(\sum^9_{n=1}\mathcal{I}^{\mu\nu\alpha\beta}_{n;(u_1,u_2,u_3,u_4)}[k_1, k_3, k_2, k_4]\Big)\Big|_{u_2 = u_4} 
\label{sumofcoefRj}
&=& - 2 \epsilon^{\mu\nu\eta}_{\ \ \ \ \kappa} (k_{1\eta}  -  k_{3\eta})
 \Big[ - G_B(u_2, u_1) + G_B(u_1, u_3) - G_B(u_3, u_2) \Big] \nonumber \\
 &\times& (k^\kappa_{2}  + k^\kappa_{4} ) \epsilon^{\alpha\beta\sigma\lambda} k_{2\sigma} k_{4\lambda}\,.
\end{eqnarray}
This expression has a very similar structure when compared to Eq. (\ref{sumcoefBj}) that we derived in the Bjorken limit.

Substituting this sum of coefficients into the expression in Eq.~(\ref{GammawithIRj}) for the box diagram in the Regge limit yields
\begin{eqnarray}
&&\Gamma^{\mu\nu\alpha\beta}_A[k_1, k_3, k_2, k_4]\Big|_{x_{Bj}\to 0} = \frac{g^2e^2 e_f^2 }{4\pi^2} \epsilon^{\mu\nu\eta}_{\ \ \ \ \kappa} (k_{1\eta}  -  k_{3\eta}) (2\pi)^4\delta^{(4)}(\sum^4_{i=1} k_i) \prod^4_{k=1}\int^1_0 du_k 
 \nonumber\\
 &&\times \Big[ - G_B(u_2, u_1) + G_B(u_1, u_3) - G_B(u_3, u_2) \Big] (k^\kappa_{2}  + k^\kappa_{4} ) \epsilon^{\alpha\beta\sigma\lambda} k_{2\sigma} k_{4\lambda} \Big[ - k_{1} \cdot k_{2} G_B(u_1, u_2)
\nonumber\\
 &&- k_{1} \cdot k_{3} G_B(u_1, u_3) - k_{1}\cdot k_{4} G_B(u_1, u_4) - k_{2} \cdot k_{3} G_B(u_2, u_3) - k_{2} \cdot k_{4} G_B(u_2, u_4) - k_{3} \cdot k_{4} G_B(u_3, u_4) \Big]^{-2}\,.
\end{eqnarray}
Using rotational invariance of the worldline, we fix $u_2 = 0$. In exact analogy to our calculation in the Bjorken limit, there are six possible orderings of the proper time variables $u_k$ and these can be split, as previously,
into the two topologies corresponding to Figs. \ref{fig:box_ordering}a and \ref{fig:box_ordering}b. As noted, the two contributions corresponding to Fig. \ref{fig:box_ordering}b have the structure $\frac{1}{(k_1+k_4)^2} \frac{1}{(k_2+k_3)^2}\rightarrow \frac{1}{(Q^2+2\,xP\cdot q)^2}$  while the four contributions corresponding to  Fig. \ref{fig:box_ordering}a are identical, having the structure $\frac{1}{l^2(Q^2+2\,xP\cdot q)}$ and contain the infrared pole in $l^2$. Since the ``cat's eye" diagrams are finite and in addition have a relative suppression $l^2/M^2$, where $M^2\equiv 2x P\cdot q\gg Q^2$ in Regge asymptotics\footnote{The longitudinal momentum fraction $x$ of the background gluons can be taken to be small but finite as $x_B\rightarrow 0$.}, we will consider only the four diagrams corresponding to Fig. \ref{fig:box_ordering}a henceforth.

Since the contributions of all for diagrams are identical, it is sufficient to examine one particular ordering of the proper time variables, which we will choose to be $u_4 > u_3 > u_1$ and multiply the expression by a factor of 4. We then get 
\begin{eqnarray}
&&\Gamma^{\mu\nu\alpha\beta}_A[k_1, k_3, k_2, k_4]\Big|_{x_{Bj}\to 0} = 4\frac{g^2e^2 e_f^2 }{4\pi^2} \epsilon^{\mu\nu\eta}_{\ \ \ \ \kappa} (k_{1\eta}  -  k_{3\eta}) (2\pi)^4\delta^{(4)}(\sum^4_{i=1} k_i) \int^1_0 du_4 \int^{u_4}_0 du_3  \int^{u_3}_0 du_1 
 \nonumber\\
 &&\times \Big[ - G_B(0, u_1) + G_B(u_1, u_3) - G_B(u_3, 0) \Big] (k^\kappa_{2}  + k^\kappa_{4} ) \epsilon^{\alpha\beta\sigma\lambda} k_{2\sigma} k_{4\lambda} \Big[ - k_{1} \cdot k_{2} G_B(u_1, 0)
\nonumber\\
 &&- k_{1} \cdot k_{3} G_B(u_1, u_3) - k_{1}\cdot k_{4} G_B(u_1, u_4) - k_{2} \cdot k_{3} G_B(0, u_3) - k_{2} \cdot k_{4} G_B(0, u_4) - k_{3} \cdot k_{4} G_B(u_3, u_4) \Big]^{-2}\,.
\end{eqnarray}

Introducing the Feynman parameters, 
\begin{eqnarray}
a_1 = 1 - u_4;\ \ \ a_2 = u_4 - u_3;\ \ \ a_3 = u_3 - u_1;\ \ \ a_4 = u_1\,,
\end{eqnarray}
and using the explicit form of the boson worldline propagator, 
\begin{eqnarray}
G_B(u_i, u_j) = |u_i - u_j| - (u_i - u_j)^2\,,
\end{eqnarray}
we can rewrite the equation as
\begin{eqnarray}
&&\Gamma^{\mu\nu\alpha\beta}_A[k_1, k_3, k_2, k_4]\Big|_{x_{Bj}\to 0} = 4\frac{g^2e^2 e_f^2 }{4\pi^2} \epsilon^{\mu\nu\eta}_{\ \ \ \ \kappa} (k_{1\eta}  -  k_{3\eta}) (2\pi)^4\delta^{(4)}(\sum^4_{i=1} k_i) \prod^4_{k=1}\int^1_0 da_k ~\delta(1 - \sum^4_{j=1}a_j)
 \nonumber\\
 &&\times \Big[ - 2a_2a_4  \Big] (k^\kappa_{2}  + k^\kappa_{4} ) \epsilon^{\alpha\beta\sigma\lambda} k_{2\sigma} k_{4\lambda} \Big[ (k_{2} + k_4)^2 a_2 a_4 + (k_{1} + k_2)^2 a_1a_3 + k^2_{1} a_3 a_4 + k^2_{2} a_1a_4 + k^2_3 a_2a_3 + k^2_4 a_1a_2 \Big]^{-2}\,.
\end{eqnarray}
where in the numerator we neglect terms\footnote{We remind the reader that due to rotational invariance of the worldline, the proper time coordinates $u_2 = 0$ and $u_2 = 1$ correspond to the same point.} proportional to $a_1 \equiv u_2 - u_4$. In  the denominator, in contrast to our discussion in the Bjorken limit, we will keep the dependence on the virtuality of the background gluons $k^2_{2}$ and $k^2_4$; these  regulate the integrals over the $\beta$ variables.

We next make the change of variables,
\begin{eqnarray}
a_1 = \alpha;\ \ \ a_i = (1-\alpha)\beta_i,\ \ \ i=2,3,4\,,
\end{eqnarray}
which yields
\begin{eqnarray}
&&\Gamma^{\mu\nu\alpha\beta}_A[k_1, k_3, k_2, k_4]\Big|_{x_{Bj}\to 0}
\\
&&= 4\frac{g^2e^2 e_f^2 }{4\pi^2} \epsilon^{\mu\nu\eta}_{\ \ \ \ \kappa} (k_{1\eta}  -  k_{3\eta}) (2\pi)^4\delta^{(4)}(\sum^4_{i=1} k_i) \int^1_0 d\alpha \prod^4_{k=2}\int^1_0 d\beta_k ~ \delta(1 - \sum^4_{j=2}\beta_j) \Big[ - 2 \beta_2 \beta_4  \Big] (k^\kappa_{2}  + k^\kappa_{4} ) \epsilon^{\alpha\beta\sigma\lambda} k_{2\sigma} k_{4\lambda}
\nonumber\\
 &&\times  \Big[ (k_{2} + k_4)^2 (1-\alpha)\beta_2 \beta_4 + (k_{1} + k_2)^2 \alpha\beta_3 + k^2_{1} (1-\alpha)\beta_3 \beta_4 + k^2_{2} \alpha \beta_4 + k^2_3 (1-\alpha)\beta_2\beta_3 + k^2_4 \alpha \beta_2 \Big]^{-2}\,.
 \nonumber
\end{eqnarray}
where in the numerator we neglect terms\footnote{These terms are of order 
$Q^2/M^2\ll 1$ in Regge asymptotics.} proportional to $\alpha = u_2 - u_4$.

The integration over $\alpha$ yields
\begin{eqnarray}
&&\Gamma^{\mu\nu\alpha\beta}_A[k_1, k_3, k_2, k_4]\Big|_{x_{Bj}\to 0} = - 2 \frac{g^2e^2 e_f^2 }{\pi^2} \epsilon^{\mu\nu\eta}_{\ \ \ \ \kappa} (k_{1\eta}  -  k_{3\eta}) (2\pi)^4\delta^{(4)}(\sum^4_{i=1} k_i) ~(k^\kappa_{2}  + k^\kappa_{4} ) \epsilon^{\alpha\beta\sigma\lambda} k_{2\sigma} k_{4\lambda}
\label{Rjfrac}\\
 &&\times  \prod^4_{k=2}\int^1_0 d\beta_k ~ \delta(1 - \sum^4_{j=2}\beta_j) \frac{\beta_2 \beta_4}{(k_{2} + k_4)^2 \beta_2 \beta_4 + k^2_{1} \beta_3 \beta_4 + k^2_3 \beta_2 \beta_3} \times \frac{1}{ (k_{1} + k_2)^2 \beta_3  + k^2_{2} \beta_4 + k^2_4 \beta_2}\,.
 \nonumber
\end{eqnarray}
We will now equate $k^2_3 = k^2_1$ and $k^2_4 = k^2_2$, which are both valid in the forward limit. As we will see, these substitutions don't affect the infrared $1/l^2$ pole in  the final answer. Hence, 
\begin{eqnarray}
&&\Gamma^{\mu\nu\alpha\beta}_A[k_1, k_3, k_2, k_4]\Big|_{x_{Bj}\to 0} = - 2 \frac{g^2e^2 e_f^2 }{\pi^2} \epsilon^{\mu\nu\eta}_{\ \ \ \ \kappa} (k_{1\eta}  -  k_{3\eta}) (2\pi)^4\delta^{(4)}(\sum^4_{i=1} k_i) ~(k^\kappa_{2}  + k^\kappa_{4} ) \epsilon^{\alpha\beta\sigma\lambda} k_{2\sigma} k_{4\lambda}
\\
 &&\times  \prod^4_{k=2}\int^1_0 d\beta_k ~ \delta(1 - \sum^4_{j=2}\beta_j) \frac{\beta_2 \beta_4}{(k_{2} + k_4)^2 \beta_2 \beta_4 + k^2_{1} \beta_3 (1 - \beta_3 )} \times \frac{1}{ (2k_1\cdot k_2 + k^2_{1})\beta_3 + k^2_2 }\,.
 \nonumber
\end{eqnarray}

The integration over the variable $\beta_3$ is dominated by small values $\beta_3\sim \frac{k^2_2}{2k_1\cdot k_2}$, where we take into account the fact that in the Regge limit $2k_1\cdot k_2 + k^2_1 \simeq 2k_1\cdot k_2$. One can therefore write,
\begin{eqnarray}
\Gamma^{\mu\nu\alpha\beta}_A[k_1, k_3, k_2, k_4]\Big|_{x_{Bj}\to 0}
 &=& - 2\, \frac{g^2e^2 e_f^2 }{\pi^2} \epsilon^{\mu\nu\eta}_{\ \ \ \ \kappa} (k_{1\eta}  -  k_{3\eta}) (2\pi)^4\delta^{(4)}(\sum^4_{i=1} k_i) ~ \frac{k^\kappa_{2}+ k^\kappa_{4}}{(k_{2} + k_4)^2 } \epsilon^{\alpha\beta\sigma\lambda} k_{2\sigma} k_{4\lambda}\nonumber \\
 &\times& \int^1_0 d\beta_3 ~ \frac{1}{ 2k_1\cdot k_2 \beta_3 + k^2_2 }\,.
\end{eqnarray}
As promised, the infrared pole $l^2 = (k_2+k_4)^2$ is manifest in this expression. 
The integration over $\beta_3$ is easily performed, giving
\begin{eqnarray}
\Gamma^{\mu\nu\alpha\beta}_A[k_1, k_3, k_2, k_4]\Big|_{x_{Bj}\to 0} 
= - 2 \frac{g^2e^2 e_f^2 }{\pi^2} \epsilon^{\mu\nu\eta}_{\ \ \ \ \kappa} (k_{1\eta}  -  k_{3\eta}) (2\pi)^4\delta^{(4)}(\sum^4_{i=1} k_i) ~ \frac{k^\kappa_{2}+ k^\kappa_{4}}{(k_{2} + k_4)^2 } \epsilon^{\alpha\beta\sigma\lambda} k_{2\sigma} k_{4\lambda}~ \frac{1}{2k_1\cdot k_2}\ln\Big[\frac{2k_1\cdot k_2 + k^2_2}{k^2_2}\Big]\,.\nonumber\\
\end{eqnarray}

Substituting this result into Eqs. (\ref{Hadgamma}) and (\ref{WModFu}), performing the integrals over intermediate coordinates and momenta, after analytical continuation to Minkowski space-time, we obtain
\begin{eqnarray}
i\tilde{W}_f^{\mu\nu}(q, P, S) = - \frac{g^2 e_f^2 }{\pi^3} \epsilon^{\mu\nu\eta}_{\ \ \ \ \kappa} q_{\eta}\,
 {\rm Im} \int \frac{d^4k}{(2\pi)^4} ~ \frac{1}{2\,k\cdot q}\ln\Big[\frac{2q\cdot k + k^2}{k^2}\Big]~  \int d^4z\, e^{-ik z} \lim_{l^\kappa\to 0} \frac{l^\kappa}{l^2 } \langle P',S| {\rm Tr_c}( F_{\alpha\beta}(z) \tilde{F}^{\alpha\beta}(0)) |P,S\rangle\,,\nonumber \\
\end{eqnarray}
where $k\equiv k_2$ is the virtuality of the background gluons, and $l = P' - P$.

We will now take into account the fact that the virtuality of the background gluons in the Regge limit is transverse, $k^2\simeq - k^2_\perp$. This is analogous to the Bjorken limit expression in Eq.~(\ref{WantisimA}), where we had $q^2 = -Q^2$. We can similarly take the imaginary part of the logarithm, which gives   
\begin{eqnarray}
&&i\tilde{W}_f^{\mu\nu}(q, P, S) =  \frac{g^2 e_f^2 }{\pi^2} \epsilon^{\mu\nu\eta}_{\ \ \ \ \kappa} q_{\eta}
 \int \frac{d^4k}{(2\pi)^4} ~ \frac{1}{2\,k\cdot q}  \int d^4z \,e^{-ik z} \lim_{l^\kappa\to 0} \frac{l^\kappa}{l^2 } \langle P',S| {\rm Tr_c}( F_{\alpha\beta}(z) \tilde{F}^{\alpha\beta}(0)) |P,S\rangle\,.
 \label{WRjbefintk}
\end{eqnarray}

The rest of the derivation follows exactly along the lines of Section~\ref{sec:Bj}.  Assuming that $k\approx (k\cdot n)P^\mu$, introducing the variable $x$ through the identity $\int dx \int d\xi e^{i(x - (k\cdot n))\xi} = 2\pi$, and performing the integration over $k$ in Eq. (\ref{WRjbefintk}), we obtain 
\begin{eqnarray}
&&i\tilde{W}_f^{\mu\nu}(q, P, S) =  \frac{4}{P\cdot q} \epsilon^{\mu\nu\eta}_{\ \ \ \ \kappa} q_{\eta} \int^1_{x_B} \frac{dx}{x} ~  \int \frac{d\xi}{2\pi} e^{-i \xi x} \lim_{l^\kappa\to 0} \frac{l^\kappa}{l^2 } \langle P',S| \frac{\alpha_s e_f^2 }{2\pi} {\rm Tr_c} F_{\alpha\beta}(\xi n ) \tilde{F}^{\alpha\beta}(0) |P,S\rangle.
\end{eqnarray}
Comparing our result with the antisymmetric part of the hadron tensor in Eq.~(\ref{WA}), we obtain finally, 
\begin{eqnarray}
S^\mu g_1(x_B, Q^2)\Big|_{x_{\rm Bj}\rightarrow 0} = \sum_f e_f^2 \frac{\alpha_s}{i\pi M_N}
 \int^1_{x_B} \frac{dx}{x} 
 \int \frac{d\xi}{2\pi} e^{-i\xi x }\lim_{l_\mu \to 0} \frac{ l^\mu }{l^2  }\langle P^\prime,S| {\rm Tr_c} F_{\alpha\beta}(\xi n) \tilde{F}^{\alpha\beta}(0) |P,S\rangle\,.
\label{resultRj}
\end{eqnarray}

\subsection{Discussion of Eq.~(\ref{resultRj})}
\label{subsec:Rj-discussion}
Eq.~(\ref{resultRj}) is the other key result of this paper. It is  valid in Regge asymptotics, where $2k\cdot q \gg Q^2 \to x > x_B$ and $x_B \to 0$. With this in mind, if we compare Eqs. (\ref{resultRj}) and (\ref{resultBj}), the result for $g_1(x_B,Q^2)$ is formally identical in the two limits. 

However both qualitatively and quantitatively, the results for $g_1(x_B,Q^2)$ can be quite different in the two limits because the scales controlling the off-forward matrix element of the topological charge density represent very different physics. In the Bjorken limit, this scale is set by $Q^2$ while in the Regge limit, it is likely set by an emergent semi-hard saturation scale $Q_s(x)$ that controls the virtuality of the background gluons $k^2 = -k_\perp^2 = -Q_s^2$. The saturation scale is a dynamical close packing scale corresponding to the maximal occupancy ($\sim 1/\alpha_s$) of gluon modes with transverse momenta $k_\perp \leq Q_s(x)$ at a given $x$~\cite{Gelis:2010nm,Kovchegov:2012mbw}. 

The scale evolution of the matrix element in Regge asymptotics will depend on the small $x$ evolution within the shock wave that describes spin diffusion from large $x$ in the polarized proton target to the DIS probe. As we noted in the introduction, there is an extensive literature studying spin diffusion in the Regge asymptotics of perturbative QCD~\cite{Kirschner:1983di,Bartels:1995iu,Bartels:1996wc,Kovchegov:2015pbl,Kovchegov:2016weo,Kovchegov:2017jxc,Chirilli:2018kkw,Boussarie:2019icw,Cougoulic:2019aja,Kovchegov:2020hgb,Cougoulic:2020tbc}.
How the physics of $U_A(1)$ breaking  discussed briefly in Sec.~\ref{subsec:Bj-discussion} is realized by small $x$ partons is an open question of great interest. Interestingly some of this nonperturbative dynamics can be explored in weak coupling since $\alpha_s(Q_s)\ll 1$ in Regge asymptotics. 
In our forthcoming work, we will write down the small $x$ effective action that is consistent with the anomalous chiral Ward identities (reflecting  the physics of the anomaly) and describe its QCD evolution. 


\section{Summary and outlook}
In this paper, we employed the worldline formalism, we developed previously for unpolarized DIS in \cite{Tarasov:2019rfp}, to compute the box diagram of polarized DIS in the Bjorken and Regge asymptotics of QCD. In particular, we find that the computation of the box diagram is nearly identical in both asymptotic limits. This is remarkable from the point of view of standard computational techniques that employ the operator product expansion, since it is well known that the OPE cannot be applied straightforwardly in Regge asymptotics~\cite{Mueller:1996hm}.  

We find that in both asymptotics, the matrix element for the  $g_1(x_B,Q^2)$ structure function is identically controlled by the triangle anomaly, which has an infrared pole in the forward scattering limit. As we discussed at length in the introduction, the cancellation of this pole involves a subtle interplay of perturbative and nonperturbative physics that is deeply related to the $U_A(1)$ problem in QCD. As a further consequence, our results bring up important questions regarding the applicability of QCD factorization to observables that are sensitive to the anomaly.  

Though the matrix element for $g_1$ is formally identical in both Bjorken and Regge limits, the underlying physics is quite different. Some of these qualitative differences are already evident from our derivations in Secs.~\ref{sec:Bj} and \ref{sec:Regge}. In paper II in this series,  we will derive the small $x$ effective action that follows from the cancellation of the infrared pole in the matrix element of the anomaly. This effective action,  consistent with anomalous chiral Ward identities, is controlled by two dimensionful scales in Regge asymptotics. The first is the color charge squared per unit area, which is proportional to the saturation scale $Q_s^2$~\cite{McLerran:1993ni,McLerran:1993ka,McLerran:1994vd}, while the second is the pure Yang-Mills topological susceptibility~\cite{Veneziano:1989ei,Shore:1990zu,Shore:1991dv}. The physics of the former has been discussed extensively in the context of small $x$ physics with the framework of the Color Glass Condensate Effective Field Theory~\cite{Gelis:2010nm}.

The latter is often tied to the physics of instantons in QCD, though this is not necessarily the case. An excellent review of the many subtleties involved can be found in~\cite{Schafer:1996wv}. More generally, one can argue that the fundamental origin of this scale has to do with the description of the  QCD vacuum as energy degenerate  $\theta$-vacua, each corresponding to distinct integer valued Chern-Simons number. We will argue that the dynamics governing helicity evolution is governed by over the barrier sphaleron transitions that are enhanced by the large dynamical saturation scale. This is analogous to the important role played by  sphaleron transitions in electroweak baryogenesis~\cite{Kuzmin:1985mm} and in QCD at  finite 
temperature~\cite{McLerran:1990de,Moore:2010jd}. Such an interplay between the saturation scale and rate of topological transitions has been studied recently for highly occupied Yang-Mills fields off-equilibrium~\cite{Mace:2016svc}.

Paper III will discuss the renormalization group evolution of the helicity dependent effective action.  With decreasing $x_B$ (or increasing energy), studies in the CGC EFT suggest that the background classical (high occupancy) fields have a lumpier 
structure with radii $\sim 1/Q_s$; as shown in \cite{Mace:2016svc} for off-equilibrium Yang-Mills fields, these lumps are also associated with a faster rate of topological transitions that flip the helicity of left handed quarks to right handed quarks or vice versa. Our conjecture, which will be explored at greater length in Paper III, is that spin diffusion in this ``disordered medium" is rapid, leading to a strong damping of $g_1(x_B,Q^2)$ and other spin-dependent observables that are sensitive to the anomaly~\cite{Shore:2007yn}. 
Because $\alpha_s(Q_s)<<1$ in Regge asymptotics, the rate of this damping can be computed in  a weak coupling framework and compared to the results of polarized DIS experiments at the EIC. Such experiments therefore have the potential to uncover novel dynamical features associated with the topology of the QCD  vacuum and its interplay with the physics of gluon saturation. 

\section*{Acknowledgements}
We are grateful to Werner Vogelsang for several very informative discussions on many aspects of the spin puzzle and for a close reading of the manuscript. We would also like to thank Gia Dvali for inspiring discussions, in particular with regard to the wider implications of our work. 
We  thank Jochen Bartels, Yoshitaka Hatta, Bob Jaffe, Jamal Jalilian-Marian, Xiangdong Ji, Yuri Kovchegov, Daniel Pitonyak, Matt Sievert, George Sterman and Feng Yuan for useful discussions on  polarized DIS relevant to this work. Fabian Rennecke and Rob Pisarski are owed thanks for their insights on instantons and the chiral anomaly. 
R.V.'s work is supported by the U.S. Department of Energy under Contract No. DE-SC0012704 and by the US DOE within the framework of the TMD Theory Topical Collaboration. His work was also supported in part by an LDRD grant from Brookhaven Science Associates. A.T.'s work is supported by the U.S. Department of Energy, Office of Science, Office of Nuclear Physics under Award No. DE-SC0004286.
