\section{Introduction}

While image-based Deep Reinforcement Learning (RL) has recently provided significant successes in various high-data domains \citep{mnih2015human,silver2017mastering,lillicrap2015continuous}, its application to physical systems remains challenging due to expensive and slow data generation, challenges with respect to safety, and the need to be robust to unexpected changes in the environment.

When training visual models in simulation, we can obtain robustness either by adaptation to target domains~\citep{ganin2016domain, bousmalis2016domain,  wulfmeier2017mutual}, or by randomising system parameters with the aim of covering all possible environment parameter changes~\citep{tobin2017domain, rajeswaran2016epopt, openai2018handmanip, pinto2017asymmetric, sadeghi2016cad2rl, viereck2017learning}. Unfortunately, training under a distribution of randomised visuals ~\citep{sadeghi2016cad2rl, viereck2017learning}, can be substantially more difficult due to the increased variability. This often leads to a compromise in final performance~\citep{openai2018handmanip, tobin2017domain}. Furthermore, it is usually not possible to cover all potential environmental variations during training. Enabling agents to generalise to unseen visuals such as distractors (task-independent aspects of the observation space) is an important question in robotics where an agent's environment is often noisy (e.g. autonomous vehicles).

To increase robustness and reduce training time, we can  make use of privileged information such as environment states, commonly accessible in simulators. By using lower-dimensional, more structured and informative representations directly as agent input, instead of noisy observations affected by visual randomisation, we can improve data efficiency and generalisation~\citep{tassa2018deepmind,peng2018sim}.

However, raw observations can be easier to obtain and dependence on privileged information during deployment can be restrictive. When exact states are available during training but not deployment, we can make use of information asymmetric actor-critic methods \citep{pinto2017asymmetric, schwab2019simultaneously} to train the critic faster via access to the state while providing only images for the actor.

By introducing Attention-Privileged Reinforcement Learning~(APRiL), we further leverage privileged information readily and commonly available during training, such as simulator states and object segmentations~\citep{todorov2012mujoco, dosovitskiy2017carla}), for increased robustness, sample efficiency, and generalisation to distractors. As a general extension to asymmetric actor-critic methods, APRiL~concurrently trains two actor-critic systems (one symmetric with a state-based agent, the other asymmetric with an image-dependent actor).
Both actors utilise attention to filter their inputs, and we encourage alignment between both attention mechanisms. As state-space learning is unaffected by visual randomisation, the observation attention module efficiently attends to state- and task-dependent aspects of the image whilst explicitly becoming invariant to task-irrelevant and noisy aspects of the environment (distractors). We demonstrate that this leads to faster image-based policy learning and increased robustness to task-irrelevant factors (both within and outside the training distribution). See Figure \ref{fig:SARL_simple_model} for a visualisation of APRiL~, its attention, and generalisation capabilities on one of our domains.

\begin{figure*}
    \centering
    \includegraphics[width=0.98\linewidth]{APRiL_simple_model.jpeg} \&
    \caption{\textbf{Model diagram (left):} APRiL~concurrently trains two attention augmented policies (one state-based, the other image-based). \textbf{Qualitative and quantitative results (middle \& right):} By aligning the observation attention to that of the state, image-based attention quickly suppresses highly varying, task-irrelevant, information \textbf{(middle second column)}. This leads to increased learning rate \textbf{(top right)} and robustness to extrapolated domains with increasing levels of unseen additional distractors \textbf{(bottom right)}. For \textit{JacoReach}, attention (\textbf{middle second column}; white and black signify high and low values) is paid only to the target object and jaco arm in training and extrapolated domains.}
    \label{fig:SARL_simple_model}
\end{figure*}

In addition, APRiL~shares a replay buffer between both agents, which further accelerates training for the image-based policy. At test-time, the image-based policy can be deployed without privileged information. We test our approach on a diverse set of simulated domains across robotic manipulation, locomotion, and navigation; and demonstrate considerable performance improvements compared to competitive baselines when evaluating on environments from the training distribution as well as in extrapolated and unseen settings with additional distractors.

\section{Problem Formulation}
\label{background}

Before introducing Attention-Privileged Reinforcement Learning~(APRiL), this section provides a
background for the RL algorithms used. For a more in-depth introduction please refer to~\citet{lillicrap2015continuous} and~\citet{pinto2017asymmetric}.

\subsection{Reinforcement Learning}

We describe an agent's environment as a Partially Observable Markov Decision Process which is represented as the tuple ${(S, O, A, P, r, \gamma, s_0)}$, where $S$ denotes a set of continuous states, $A$ denotes a set of either discrete or continuous actions, $P: S \times A \times S \rightarrow \{ x \in \mathbb{R} | 0 \le x \le 1\}$ is the transition probability function, $r : S \times A \rightarrow \mathbb{R}$ is the reward function, $\gamma$ is the discount factor, and $s_0$ is the initial state distribution. $O$ is a set of continuous observations corresponding to continuous states in $S$. At every time-step $t$, the agent takes action $a_t = \pi(\cdot|s_t)$ according to its policy $\pi : S \rightarrow A$. The policy is optimised as to maximize the expected return $R_t = E_{s_0}[\sum_{i=t}^{\infty}\gamma^{i-t}r_i|s_0]$. The agent's Q-function is defined as $Q_\pi(s_t, a_t) = E[R_t|s_t, a_t]$.

\subsection{Asymmetric Deep Deterministic Policy Gradients}
\label{sec:addpg}

\begin{wrapfigure}{r}{0.45\textwidth}
    \centering
    \vspace{-20pt}
    \includegraphics[width=0.98\linewidth]{APRiL_diagram_2.jpeg}
    \caption{\textbf{APRiL~'s architecture}. Blue, green and orange represent symmetric and asymmetric actor critic and attention alignment modules ($A_s$, $A_o$, $A_T$). The diamond represents the attention alignment loss. Dashed and solid blocks are non-trainable and trainable networks. The $\otimes$ operator signifies element-wise multiplication. Experiences are shared using a shared replay buffer.}
    \vspace{-60pt}
    \label{fig:SARL_MODEL}
\end{wrapfigure}

Asymmetric Deep Deterministic Policy Gradients (asymmetric DDPG)~\citep{pinto2017asymmetric} represents a type of actor-critic algorithm designed specifically for efficient learning of a deterministic, observation-based policy in simulation. This is achieved by leveraging access to more compressed, informative environment states, available in simulation, to speed up and stabilise training of the critic.

The algorithm maintains two neural networks: an observation-based actor or policy $\pi_{\theta} : O \rightarrow A$ (with parameters ${\theta}$) used during training and test time, and a state-based Q-function (also known as critic) $Q^\pi_{\phi} : S \times A \rightarrow R$ (with parameters ${\phi}$) which is only used during training.

To enable exploration, the method (like its symmetric version~\citep{silver2014deterministic}) relies on a noisy version of the policy (called behavioural policy), e.g. $\pi_b(o) = \pi(o) + z$ where $z \sim \mathcal{N}(0, 1)$ (see Appendix \ref{hyperparams} for our particular instantiation). The transition tuples $(s_t, o_t, a_t, r_t, s_{t+1}, o_{t+1})$ encountered during training are stored in a replay buffer~\citep{mnih2015human}. Training examples sampled from the replay buffer are used to optimize the critic and actor. By minimizing the Bellman error loss $\mathcal{L}_{critic} = (Q(s_t, a_t) - y_t)^2$, where $y_t = r_t+\gamma Q(s_{t+1}, \pi(o_{t+1}))$, the critic is optimized to approximate the true Q values. The actor is optimized by minimizing the loss: \newline $\mathcal{L}_{actor} = - E_{s, o \sim \pi_b(o)}[Q(s, \pi(o))]$.

\section{Attention-Privileged Reinforcement Learning~(APRiL)}

APRiL~improves the robustness and sample efficiency of an observation-based agent by using multiple ways to benefit from privileged information. First, we use an asymmetric actor-critic setup~\citep{pinto2017asymmetric} to train the observation based actor. Second, we additionally train a quicker learning state-based actor, while sharing replay buffers, and aligning attention mechanisms between both actors. We emphasise here that our approach can be applied to any asymmetric, off-policy, actor-critic method~\citep{konda2000actor} with the expectation of similar performance benefits to those demonstrated in this paper. Specifically we choose to build off Asymmetric DDPG~\citep{pinto2017asymmetric} due to its accessibility.

APRiL~is comprised of three modules as displayed in Figure \ref{fig:SARL_MODEL}. The first two modules, ${A_s}$ and ${A_o}$, each represent a separate actor-critic with an attention network incorporated over the input for each actor. For the \textit{state-based} module ${A_s}$ we use standard symmetric DDPG, while the \textit{observation-based} module ${A_o}$ builds on asymmetric DDPG, with the critic having access to states. Finally, the third part ${A_\target}$ represents the alignment process between attention mechanisms of both actor-critic agents to more effectively transfer knowledge between the both learners respectively.

\textbf{${A_s}$} consists of three networks: $Q^\pi_s$, $\pi_s$, ${h}_s$ (critic, actor, and attention) with parameters $\{{{\phi}_s}, {{\theta}_s}, {{\psi}_s}\}$. Given input state $s_t$, the attention network outputs a soft gating mask ${h}_t$ of same dimensionality as the input, with values ranging between $[0,1]$. The input to the actor is an attention-filtered version of the state, $s^a_t = {h}_s(s_t) \odot s_t$. To encourage a sparse masking function, we found that training this attention module on both the traditional DDPG loss as well as an entropy loss helped:
\begin{equation}
  \mathcal{L}_{h_s} = - E_{s \sim \pi_b}[Q_s(s, \pi_s(s^a)) - \beta H({h}_s(s))],
  \label{eq: actor_state_loss}
\end{equation}
    \noindent where $\beta$ is a hyperparameter (set through grid-search, see Appendix \ref{hyperparams}) to weigh the additional entropy objective, and $\pi_b$ is the behaviour policy that obtained experience (in this case from a shared replay buffer). The actor and critic networks $\pi_s$ and $Q_s$ are trained with the symmetric DDPG actor and Bellman error losses. We found that APRiL~was not sensitive to the absolute value of $\beta$, only the magnitude, and was set low enough to not suppress task-relevant parts of the state-space.

Within \textbf{${A_\target}$}, the state-attention obtained in ${A_s}$ is converted to corresponding observation-attention ${T}$ to act as a self-supervised target for the observation attention module in ${A_o}$.
This is achieved in a two-step process. First, state-attention ${h}_s(s)$ is converted into object-attention $c$, which specifies how task-relevant each object in the scene is. The procedure uses information about which dimension of the environment state relates to which object. Second, object-attention is converted to observation-space attention by performing a weighted sum over object-specific segmentation maps \footnote{Simulators (e.g., \citep{todorov2012mujoco, dosovitskiy2017carla}) commonly provide functionality to access these segmentations and semantic information for the environment state.}:
\begin{equation}
    {c} = M \cdot {h}_s(s), \qquad {T} =  \sum_{i=0}^{N - 1} {c}_i \cdot {z}_i
    ~~~~\label{eqn_3}
\end{equation}

Here, $M \in \{0,1\}^{N \times n_s}$ ($n_s$ is the dimensionality of $s$) is an environment-specific, predefined adjacency matrix that maps the dimensions of $s$ to each corresponding object, and ${c} \in [0,1]^N$ is an attention vector over the $N$ objects in the environment. ${c}_i$ corresponds to the $i^{th}$ object attention value. ${z}_i \in \{0,1\}^{W \times H}$ is the binary segmentation map of the $i^{th}$ object segmenting the object with the rest of the scene, and has the same dimensions as the image. ${z}_i$ assigns values of $1$ for pixels in the image occupied by the $i^{th}$ object, and $0$ elsewhere. $T \in [0,1]^{W \times H}$ is the converted state-attention to observation-space attention to act as a target on which to train the observation-attention network ${h}_o$.
\begin{algorithm}[h!]
  \caption{Attention-Privileged Reinforcement Learning}
  \label{alg:example}
\begin{algorithmic}
  \STATE Initialize the actor-critic modules ${A_s}$, ${A_o}$, attention alignment module $A_T$, replay buffer $R$
  \FOR{episode$=1$ {\bfseries to} $M$}
  \STATE Initial state $s_0$
  \STATE Set DONE $\leftarrow$ FALSE
  \WHILE{$\neg$ DONE}
  \STATE Render image observation $o_t$ and segmentation maps $z_t$:
  \STATE \hskip4.5em  $o_t, z_t \leftarrow$ {renderer}$(s_t)$
  \IF{episode \textbf{mod} $2 = 0$}
  \STATE Obtain action $a_t$ using obs-behavioral policy and obs-attention network:
  \STATE \hskip4.5em  $a_t \leftarrow \pi_{o}(h_o(o_t) \odot o_t)$
  \ELSE
  \STATE Obtain action $a_t$ using state-behavioral policy and state-attention network:
  \STATE \hskip4.5em  $a_t \leftarrow \pi_{s}(h_s(s_t) \odot s_t)$
  \ENDIF
  \STATE Execute action $a_t$, receive reward $r_t$, DONE flag, and transition to $s_{t+1}$
  \STATE Store $(s_t, o_t, z_t, a_t, r_t, s_{t+1}, o_{t+1})$ in $R$
  \ENDWHILE
  \FOR{$n=1$ {\bfseries to} $N$}
  \STATE Sample minibatch $\{s, o, z, a, r, s^{'}, o^{'}\}_0^B$ from $R$
  \STATE Optimise state- critic, actor, and attention using $\{s, a, r, s^{'}\}_0^B$ with ${A_s}$
  \STATE Convert state-attention to target observation-attention $\{T\}_0^B$ using $\{s, o, z\}_0^B$ with $A_T$
  \STATE Optimise observation- critic, actor, and attention using $\{s, o, T, a, r, s^{'}, o^{'}\}_0^B$ with ${A_o}$
  \ENDFOR
  \ENDFOR
\end{algorithmic}
\label{algo:APRiL}
\end{algorithm}
The observation module ${A_o}$ also consists of three networks: $Q^\pi_o$, $\pi_o$, ${h}_o$ (respectively critic, actor, and attention) with parameters $\{{\phi_o},{{\theta}_o}, {{\psi}_o}\}$. The structure of this module is the same as ${A_s}$ except the actor and critic now have asymmetric inputs. The actor's input is the attention-filtered version of the observation, $o^a_t = {h}_o(o_t) \odot o_t$ \footnote{In practice, the output of ${h}_o(o_t)$ is tiled to match the number of channels that the image contains}. The actor and critic $\pi_o$ and $Q_o$ are trained with the asymmetric DDPG actor and Bellman error losses in \ref{sec:addpg}.
The main difference between ${A_o}$ and ${A_s}$ is that the observation attention network ${h}_o$ is trained on both the actor loss and an object-weighted mean squared error loss:
\begin{equation}
\mathcal{L}_{h_o} = E_{o, s \sim \pi_b} [\frac{1}{2}\sum_{ij}\frac{1}{w_{ij}}({h}_o(o) - T)^2_{ij}  - \nu Q_o(s, \pi_o(o^a))]
\end{equation}

where weights $w_{ij}$ denote the fraction of the image $o$ that the object present in $o_{i,j,1:3}$ occupies, and $\nu$ represents a hyperparameter for the relative weighting of both loss components (see Appendix \ref{hyperparams} for exact value). The weight terms, $w$, ensure that the attention network becomes invariant to the size of objects during training and does not simply fit to the most predominant object in the scene.

During training, experiences are collected evenly from both state and observation based agents and stored in a shared replay buffer (similar to~\citet{schwab2019simultaneously}). This is to ensure that: 1. Both state-based critic $Q_s$ and observation-based critic $Q_o$ observe states that would be visited by either of their respective policies. 2. The attention modules $h_s$ and $h_o$ are trained on the same data distribution to better facilitate alignment. 3. Efficient discovery of highly performing states from $\pi_s$ are used to speed up learning of $\pi_o$.

Algorithm \ref{algo:APRiL} shows the pseudocode for a single actor implementation of APRiL. In practice, in order to speed up data collection and gradient computation, we parallelise the agents and environments and ensure equal data is generated by state- and image-based agents.

\section{Experiments}

\begin{figure*}
    \includegraphics[width=0.32\linewidth]{NavIncDDPGv2.png}
    \includegraphics[width=0.32\linewidth]{JacIncDDPG.png}
    \includegraphics[width=0.32\linewidth]{WalkIncDDPG.png} \\
    \includegraphics[width=0.32\linewidth]{NavAblNewName.png}
    \includegraphics[width=0.32\linewidth]{JacAbl2Large.png}
    \includegraphics[width=0.32\linewidth]{WalkAblLarge.png}

    \caption{Learning curves for observation-based policies during Domain Randomisation (DR). \textbf{Top row}: comparison with baselines. \textbf{Bottom row}: comparison with ablations. \textbf{Solid line}: mean performance. \textbf{Shaded region}: covers minimum and maximum performances across $5$ seeds. APRiL's attention and shared replay lead to stronger or commensurate performance.}
    \label{fig:SARL_experiments}
\end{figure*}

\begin{figure*}
    \centering
    \includegraphics[width=0.36\linewidth]{NavHeldOutDDPG.png}
    \includegraphics[width=0.36\linewidth]{JacHeldOutDDPG.png}
    \includegraphics[width=0.18\linewidth]{WalkHeldOutDDPG.png}
    \caption{Comparing average return of the image-policy between training, interpolated and extrapolated domains (100 each). Plots reflect mean and 2 standard deviations for average return (5 seeds). APRiL~generalises due to its attention and outperforms the baselines. We compare against a random agent to gauge the degree of degradation in policy performance between domains.}
    \label{fig:percent_decrease}
\end{figure*}
\begin{figure*}
    \centering
    \includegraphics[width=0.99\linewidth]{SARL_video_sub_large_2001.jpeg}
    \caption{Held-out domains and APRiL~attention maps. For the extrapolated domain columns (extra), top and bottom represent \textbf{ext-4} and \textbf{ext-8}. White/black signify high/low attention values. Attention suppresses the background and distractors and helps generalise.}
    \label{fig:hold_out_envs}
\end{figure*}

\label{sec:protocol}

We evaluate APRiL~over the following environments (see Appendix~\ref{app:environments} for more details): {1. \textit{NavWorld}: the circular agent is sparsely rewarded for reaching the triangular target in the presence of distractors.} {2. \textit{JacoReach}: the Kinova arm is rewarded for reaching the diamond-shaped object in the presence of distractors.} {3. \textit{Walker2D}: this slightly modified (see Appendix~\ref{app:environments}) Deepmind Control Suite environment~\citep{tassa2018deepmind} the agent is rewarded for walking forward whilst keeping its torso upright.}

During training, for APRiL~, its ablations, and all baselines, we perform Domain Randomisation (DR) \cite{tobin2017domain, sadeghi2016cad2rl}, randomising the following environment parameters to enable generalisation with respect to them: camera position, orientation, textures, materials, colours, object locations, background (see Appendix~\ref{app:rand}).

We start by comparing APRiL~against two competitive baselines that also exploit privileged information during training. We compare against the \textit{Asymmetric DDPG} (asym-DDPG) baseline \citep{pinto2017asymmetric} to evaluate the importance of privileged attention and shared replay for learning and robustness to distractors. Our second baseline, \textit{State-Mapping Asymmetric DDPG} (s-map asym-DDPG), introduces a bottleneck layer trained to predict the environment state using an \textit{$L_{2}$} loss. This is another intuitive approach that further exploits state information in simulation~\citep{zhang2016vision} to learn informative representations that are robust to visual randomisation. This approach does not incorporate object-centric attention or leverage privileged object segmentations. We note that since this baseline learns state estimation it is not expected to extrapolate well to domains with additional distractor objects and varying state spaces (with respect to the training domain). We also compare APRiL~with DDPG to emphasise the difficulty of these DR tasks if privilege information is not leveraged.

We perform an ablation study to investigate which components of APRiL~contribute to performance gains. The ablations consist of: {1. \textit{APRiL~no sup}: the full setup except without attention alignment. Here the observation attention module must learn without guidance from the state-agent.} {2. \textit{APRiL~no share}: APRiL~without a shared replay.} {3. \textit{APRiL~no back}: uniform object attention values $c$ are used to train the observation attention module, thereby only suppressing the background. Here we investigate the importance of object-suppression for generalisation.}

We investigate the following to evaluate how well APRiL~facilitates transfer across visually distinct domains: Does~APRiL~:
{1. Increase \textbf{sample-efficiency} during training?}
{2. Affect \textbf{interpolation} performance on unseen environments from the training distribution?}
{3. Affect \textbf{extrapolation} performance on environments outside the training distribution?}

\subsection{Performance On The Training Distribution}

Figure \ref{fig:SARL_experiments} shows that APRiL~outperforms the baselines for each environment (except \textit{Walker2D} where it matches s-map asym-DDPG). The ablations in Figure \ref{fig:SARL_experiments} show that a shared replay buffer, background suppression, and attention-alignment each individually provide benefits but are most effective when combined together. Interestingly, background suppression is extremely effective for sample-efficiency, as for these domains the majority of the irrelevant, highly varying, aspects of the observation-space are occupied by the background. It is also surprising that the s-map asym-DDPG baseline, which learns to map to environment states, does not outperform asym-DDPG and does not match APRiL~'s
performance for \textit{NavWorld} and \textit{JacoReach}. For these domains, predicting states (including those of distractors) is difficult\footnote{For \textit{JacoReach} prediction errors and policy performance are sensitive to state-space. In Figure \ref{fig:SARL_experiments} we plot the best performing state-space. Refer to Appendix \ref{app:s-map-abl} for further details.} and prediction errors limit policy performance. For \textit{Walker2D}, in the absence of distractor objects, s-map asym-DDPG is a competitive baseline and APRiL~provides
marginal gains.

\subsection{Interpolation: Transfer To Domains From The Training Distribution}
\label{interp}

We evaluate performance on environments unseen during training but within the training distribution (see Appendix \ref{app:rand}). For \textit{NavWorld} and \textit{JacoReach}, the interpolated environments have the same number of distractors, sampled from the same object catalogue, as the training distribution. Figure \ref{fig:percent_decrease} plots the return on these held-out domains. For all algorithms, we observe minimal degradation in performance between training and interpolated domains. However, as APRiL~outperforms on the training distribution (apart from; \textit{Walker2D} for s-map asym-DDPG, \textit{JacoReach} for APRiL~no back), its final performance on the interpolated domains is significantly better, emphasising the benefits of both privileged attention and a shared replay.

\subsection{Extrapolation: Transfer To Domains Outside The Training Distribution}
\label{extrap}

For \textit{NavWorld} and \textit{JacoReach}, we investigate how well each method generalises to extrapolated domains with additional distractor objects (specifically 4 or 8; referred as \textbf{ext-4} and \textbf{ext-8}). The textures and colours of these objects are sampled from a held-old out set not seen during training. The locations are sampled; randomly for \textit{NavWorld}, from extrapolated arcs of two concentric circles of different radii for \textit{JacoReach}. Shapes are sampled from the training catalogue of distractors. We do not extrapolate for \textit{Walker2D}, as this domain does not contain distractors. However, we show (in the previous sections) that APRiL~is still beneficial during DR for this domain and therefore demonstrate its application does not need to be restricted to environments with clutter. Please refer to Figure \ref{fig:hold_out_envs} for examples of the extrapolated domains.

Figure \ref{fig:percent_decrease} shows that APRiL~generalises and performs considerably better on the held-out domains than each baseline. Specifically, when comparing with the baselines that leverage privilege information, for \textit{JacoReach} performance falls by $11$\%\footnote{Percentage decrease is taken with respect to additional return over a random agent
on the training domain.} for APRiL~instead of $42$\% and $48$\% for asym-DDPG and s-map asym-DDPG respectively. The ablations demonstrate that effective distractor suppression is crucial for generalisation. This is particularly prominent for \textit{JacoReach} where the performance drop for the methods that use attention-alignment (APRiL~and APRiL~no share) is $11$\% and $15$\%, which is far less than $27$\% and $51$\% (APRiL~no back and APRiL~no sup) for those that do not learn to effectively suppress distractors.

\subsection{Attention Module Analysis}
\label{attention_analysis}

We visualise APRiL's attention maps (Figure \ref{fig:hold_out_envs}, \ref{fig:ATTMAPS_2}, \ref{fig:ATTMAPS} (in Appendix \ref{app:s-map-abl}) and these \href{https://sites.google.com/view/april-domain-randomisation/home}{videos}) on both interpolated and extrapolated domains. For \textit{NavWorld}, attention is correctly paid to all relevant aspects (agent and target; circle and triangle respectively) and generalises well. For \textit{JacoReach}, attention suppresses the distractors even on the extrapolated domains, achieving robustness with respect to them. Interestingly, as we encourage sparse attention, APRiL~learns to only pay attention to every-other-link of the arm (as the state of an unobserved link can be inferred by observing those of the adjacent links). For \textit{Walker2D}, dynamic object attention is learnt (different objects are attended based on the state of the system - see Figure \ref{fig:ATTMAPS}). When upright, walking, and collapsing, APRiL~pays attention to the lower limbs, every other link, and foot and upper body, respectively. We suspect that in these scenarios, the optimal action depends most on the state of the lower links (due to stability), every link (coordination), and foot and upper body (large torque required), respectively.

\section{Related Work}
A large body of work investigates the problem of learning robust policies that generalise well outside of the training distribution.
Work on \textbf{transfer learning} leverages representations from one domain to efficiently solve a problem from a different domain~\citep{donahue2014decaf, oquab2014learning, rusu2016sim}.
In particular, \textbf{domain adaptation} techniques aim to adapt a learned model to a specific target domain, often optimising models such that representations are invariant to the shift in the target domain ~\citep{ganin2016domain, long2015learning, bousmalis2016domain, wulfmeier2017mutual}.
These methods commonly require data from the target domain in order to transfer and adapt effectively.

In contrast, \textbf{domain randomisation} covers a distribution of environments by randomising visual~\citep{tobin2017domain} or dynamical parameters~\citep{peng2018sim} during training in order to generalise~\citep{sadeghi2016cad2rl, rajeswaran2016epopt, viereck2017learning,held2017probabilistically, openai2018handmanip, laskin2020reinforcement}.
In doing so, such methods shift the focus from adaptation to specific environments to \textbf{generalisation and robustness} by covering a wide range of variations.
Recent work automatically varies this distribution during training~\citep{akkaya2019solving} or trains a canonical invariant image representation~\citep{james2019sim}.
However, while randomisation can enable us to learn robust policies, it significantly increases training time due to the increased environment variability~\citep{openai2018handmanip}, and can reduce asymptotic performance. Our work partially addresses this fact by training two agents, one of which is not affected by visual randomisations.

Other works explicitly encourage representations \textbf{invariant} to observation space variations \citep{slaoui2019robust, srinivas2020curl, james2019sim}. Contrastive techniques \citep{slaoui2019robust, srinivas2020curl} use a clear separation between positive (and negative) examples, predefined by the engineer, to encourage \textbf{invariance}. Unlike APRiL, these \textbf{invariances} are over abstract spaces and are not designed to exploit \textbf{privileged information}; shown to be beneficial by APRiL's ablations. Furthermore, APRiL's \textbf{invariance} is \textbf{task-driven} via attention. Approaches like \citep{james2019sim,zhang2016vision}, learn invariances in a supervised manner, mapping from observation to a predefined space. Unlike APRiL, these methods are unable to discover task-independent aspects of the mapping-space, limiting robustness and generalisation. Finally, unlike APRiL~as a model-free RL approach, some model-based works use forward or inverse models \citep{pathak2017curiosity, hafner2019learning, schmeckpeper2019learning} to achieve invariance.

Existing comparisons in the literature demonstrate that, even without domain randomisation, the increased dimensionality and potential partial observability complicates learning for RL agents~\citep{tassa2018deepmind,schwab2019simultaneously}. In this context, accelerated training has also been achieved by using \textbf{access to privileged information} such as environment states to asymmetrically train the critic in actor-critic RL~\citep{schwab2019simultaneously,pinto2017asymmetric,foerster2018counterfactual}.
In addition to using additional information to train the critic, \citep{schwab2019simultaneously} use a \textbf{shared replay buffer} for data generated by image- and state-based actors to further accelerate training for the image-based agent.
Our method extends these approaches by sharing information about relevant objects by aligning agent-integrated attention mechanisms between an image- and state-based actors.

Recent experiments have demonstrated the strong dependency and interaction between attention and learning in human subjects~\citep{LEONG2017451}. In the context of machine learning, \textbf{attention mechanisms} have been integrated into RL agents to increase robustness and enable interpretability of an agent's behaviour~\citep{sorokin2015deep, Mott2019TowardsIR}. In comparison, we focus on utilising the attention mechanism as an interface to transfer information between two agents to enable faster training and better generalisation.

\section{Conclusion}
\label{conclusion}

We introduce Attention-Privileged Reinforcement Learning~(APRiL), an extension to asymmetric actor-critic algorithms that leverages attention mechanisms and access to privileged information such as simulator environment states. The method benefits in two ways in addition to asymmetry between actor and critic: via aligning attention masks between image- and state-space agents, and by sharing a replay buffer.
Since environment states are not affected by visual randomisation, we are able to learn efficiently in the image domain especially during domain randomisation where feature learning becomes increasingly difficult.
Evaluation on a diverse set of environments demonstrates significant improvements over competitive baselines including asym-DDPG and s-map asym-DDPG; and show that APRiL~learns to generalise favourably to environments not seen during training (both within and outside of the training distribution). Finally, we investigate the relative importance of the different components of APRiL~ in an extensive ablation.


\section{Introduction}

Over the last decades, our understanding of hadron structure has remarkably improved
thanks to impressive experimental and theoretical progress.
In particular,  there is currently much  activity aiming at further unravelling
the nucleon's spin structure.

As is well known, the total quark and
anti-quark  spin contribution to the nucleon spin
was found to be only about $\sim 25 \%$\cite{deFlorian:2014yva,Nocera:2014gqa,Sato:2016tuz},  information that has been mostly provided
by  Deep-Inelastic scattering (DIS) data.

One emphasis is on the determination of the
spin-dependent gluon distribution, $\Delta g$,
of the nucleon, which ultimately would give the gluon  contribution to
the nucleon spin.
Clear evidence for a non-vanishing polarization of gluons was found in the region of momentum fraction and at the scales mostly probed by the BNL Relativistic Heavy Ion Collider (RHIC) data \cite{deFlorian:2014yva}. But from the available results it is not possible yet to set a precise value on the total gluon contribution to the proton spin and  to the flavor decomposition of the quarks \footnote{See \cite{deFlorian:2019egz} for an analysis of the proton spin budget at three loops.}.

Spin asymmetries in
high-energy $pp$ scattering can be particularly sensitive to $\Delta g$,
for processes where gluons in the initial state contribute already at
the lowest order of perturbation theory. Furthermore, the Drell-Yan process, the production
of vector bosons \cite{Adare:2018csm,Adam:2018bam}, has been shown to be very relevant to help for quark flavor separation,
 since quarks couple differently to both $W^{\pm}$ and $Z$ bosons \cite{deFlorian:2010aa}.

Therefore, one particularly interesting process that satisfies both requirements
is the single-inclusive production of large transverse-momentum
($q_T$) vector bosons, $pp \to V X$, that can help to provide information on both $\Delta g$ and the flavour decomposition of the quark contribution.

In order to make reliable quantitative predictions for a high-energy process, it is crucial to determine the next-to-leading order (NLO) QCD corrections to the Born approximation. In general, the key issue here is to check the perturbative stability of the process considered, i.e. to examine to what extent  the NLO corrections affect the cross sections and  spin asymmetries relevant for experimental measurements. Only if the corrections are under control can a process that shows good sensitivity to, say, $\Delta g$ at the lowest order be regarded as a genuine probe of the polarized gluon distribution and be reliably used to extract it from future data.

In the particular case of (large transverse momentum) gauge boson production another  perturbative issues appear due to the existence of two physical scales, the transverse momentum $q_T$ and the `mass' $Q$ of the boson (or its virtuality, specially in the case of  photon production).

In the large-$q_T$ region ($q_T\sim Q$), where the transverse momentum is of the order of the vector boson mass, the QCD perturbative series is controlled by a small expansion parameter, $\alpha_s(Q)$, and  calculations based on the truncation of the perturbative series at a fixed order in $\alpha_s$
are theoretically justified. In this region, the QCD radiative corrections for the unpolarized cross section are known up to the next-to-leading order (NLO) \cite{Ellis:1981hk,Arnold:1988dp,Gonsalves:1989ar} in an  analytical form and next-to-next-to-leading order (NNLO) corrections were recently obtained in numerical implementations in \cite{Ridder:2016nkl}--\cite{Gehrmann-DeRidder:2017mvr}. In the polarized case, a first attempt to achieve the cross section at NLO accuracy for virtual photon production was presented in \cite{Chang:1997ik}, where only the non-singlet contribution was obtained. The full NLO result for the double polarized cross section, including the relevant gluon initiated channels, was missing so far.

Nonetheless the bulk of the vector boson events (particularly at RHIC with a center-of-mass energy of $\sqrt{S}=510$ GeV) is produced in the small-$q_T$ region ($q_T\ll Q$), where the convergence of the fixed-order expansion is spoiled by the presence of large logarithmic terms, $\alpha_s^n\ln^m (Q^2/q_T^2)$. To obtain reliable predictions, these logarithmically-enhanced terms
have to be systematically resummed to all perturbative orders (see e.g.\cite{Catani:2015vma},\cite{Bizon:2019zgf} and references therein). But even in that case, the fixed-order calculation becomes essential in order to perform the proper matching with the resummed contribution, affecting the full result even at rather small transverse momentum. Furthermore, since the small-$q_T$ region is mostly affected by soft gluon emission, which is independent on the polarization of the emitting parton, the resummation typically affects both polarized and unpolarized cross section in a rather similar way, rendering a very small effect at the level of asymmetries.

Therefore, counting with the NLO corrections for polarized hadronic collisions becomes fundamental in both kinematical regimes in order to understand the data produced at RHIC and to extract the corresponding information in terms of polarized partonic distributions. In this paper we compute the NLO corrections to the cross section for the production of gauge bosons ($Z$, $W^\pm$ and $\gamma^*$) with large transverse momentum in double (longitudinally) polarized hadronic collisions.

The paper is organized as follows. In Sect.~\ref{sect1} we briefly review the procedure to compute the NLO corrections with emphasis on the issues arising due to the treatment of $\gamma^5$ in dimensional regularization. In Sect.~\ref{sect2} we present the analytical results for the NLO corrections for gauge boson  production (some of the lengthly functions are given in the Appendix). In Sect.~\ref{sect3} we study the phenomenological impact of the NLO corrections at RHIC kinematics, including the measurable double longitudinal asymmetries. Finally, in Sect.~\ref{sec:summa} we summarize our results.

\section{Cross Section Calculation}
\label{sect1}

\subsection{Inclusive cross section}

We calculate the inclusive cross section of electroweak boson production in the framework of perturbative QCD, considering the case of both polarized and unpolarized initial hadrons. The process is described as

\begin{equation}
 h_1(P_1)+h_2(P_2) \rightarrow V(Q)+X,
\end{equation}

\noindent where $h_i$, $i=1,2$ are the polarized/unpolarized initial hadrons with momenta $P_i$, and $V$ is either a $W^{\pm}$, $Z^0$ or virtual photon $\gamma^*$ with $Q^2=M_V^2$, energy $E_Q=Q^0$ and large transverse momentum $q_T$ with respect to the collision axis. The unpolarized $\sigma$ and polarized $\Delta \sigma$ cross sections are obtained from the sum or difference of helicity-dependent cross sections given by

\begin{equation}
 \sigma = \frac{1}{4} \left( \sigma^{++} + \sigma^{+-} + \sigma^{-+} + \sigma^{--} \right),
\end{equation}

\noindent for the unpolarized case, and

\begin{equation}
 \Delta \sigma = \frac{1}{4} \left( \sigma^{++} - \sigma^{+-} - \sigma^{-+} + \sigma^{--} \right),
\end{equation}

\noindent for the polarized one. Here the superindices $+,-$ denote the helicities of the two incoming hadrons.

The inclusive (polarized) unpolarized cross section can be expressed as

\begin{equation}
 E_Q \frac{d (\Delta) \sigma}{d^3 Q} = \sum_{a,b} \int^1_0 dx_a \ dx_b \ (\Delta)f_a^{h_1}(x_a,\mu_F^2) \ (\Delta)f_b^{h_2}(x_b,\mu_F^2) \ E_Q \frac{d (\Delta) \hat{\sigma}^{a,b}}{d^3 Q}(p_a, p_b, \mu_F^2).
\end{equation}

\noindent Here $(\Delta)f^h_a(x, \mu_F^2)$ is the (polarized) unpolarized parton distribution function (PDF) of parton $a$ with momentum fraction $x$ in hadron $h$, probed at the scale $\mu_F^2$, which are given by the combinations

\begin{equation}
 f_a(x,\mu_F^2)=f^{+}_a(x,\mu_F^2)+f_a^{-}(x,\mu_F^2), \quad \Delta f_a(x,\mu_F^2)=f_a^{+}(x,\mu_F^2)-f_a^{-}(x,\mu_F^2).
\end{equation}

\noindent In this case, the superindex $+,-$ denotes de helicity orientation of the parton $a$ with respect to the helicity of the parent hadron. The perturbative (polarized) unpolarized cross section $(\Delta) \hat{\sigma}^{a,b}(p_a, p_b, \mu_F^2)$ corresponds to the hard-scattering partonic process

\begin{equation}
 a(p_a)+b(p_b) \rightarrow V(Q)+X,
\end{equation}

\noindent where $p_i$ represents the parton momentum. Here collinear singularities due to the radiation of massless partons are factorized out at the scale $\mu_F^2$ and included in the scale dependent (polarized) unpolarized PDFs $(\Delta)f^h_a(x, \mu_F^2)$.

The partonic momenta can be described in terms of the hadronic ones using the relation $p_i=x_i P_i$. As it is customary, we introduce the Mandelstam variables for both the hadronic and partonic levels
\begin{equation}
    \begin{aligned}[c]
        S \equiv &(P_1+P_2)^2, \quad T \equiv (P_1-Q)^2, \quad U \equiv (P_2-Q)^2, \\
        s \equiv &(p_a+p_b)^2, \quad t \equiv (p_a-Q)^2, \quad u \equiv (p_b-Q)^2, \\
        S_{23}& \equiv S+T+U-Q^2, \quad s_{23} \equiv s+t+u-Q^2.
    \end{aligned}
\end{equation}

\noindent $S$ and $s$ are the hadronic and partonic invariant center-of-mass energies squared of the colliding system, respectively, while $s_{23}$ is the invariant mass squared of the system recoiling against the boson $V$. The partonic cross section $\hat{\sigma}^{a,b}$ has a singular behaviour in the limit $s_{23} \rightarrow 0$ related to the cancellation of singularities due to soft gluon emission in the recoiling system and virtual gluon infrared singularities. However, the cross section is actually integrable and the  $1/s_{23}$ terms can be dealt with the following change of variables in the momentum fraction integration

\begin{equation}
    \int^1_0 \int^1_0 dx_1 \ dx_2 \ \theta(s_{23}) \ \theta(p_1^0+p_2^0-E_Q)=\int^1_B dx_1 \ \frac{1}{x_1 \ S + U - Q^2} \ \int^A_0 ds_{23},
\label{eq_ints23}
\end{equation}

\noindent where the new integration limits are given by
\begin{equation}
    A=U+x_1 \ (S_{23}-U), \quad B=-\frac{U}{S_{23}+U}.
\end{equation}

\subsection{Perturbative calculation techniques}

The analytical procedures to compute the perturbative cross section are already well established. At the lowest order (LO) only two channels contribute, as depicted in Fig.\ref{fig_FeynBorn}: the annihilation process $q \overline{q} \rightarrow V g$ (a) and the Compton process $q g \rightarrow V q$ (b). At the following order (NLO) virtual and real contributions must be considered. In this case, there are four contributing channels: the process initiated by a quark and an antiquark $q \overline{q} \rightarrow V (q \overline{q}, \ q' \overline{q}', \ gg)$, the one initiated by two quarks $q q \rightarrow V q q$, the process initiated by a single gluon $q g \rightarrow V q g$ and the two gluon process $g g \rightarrow V q \overline{q}$.

\begin{figure}
\centering
  \includegraphics[width=0.4\linewidth]{Fig_Feynman_Born.pdf}
  \caption{Born diagrams of the two processes at LO.}
  \label{fig_FeynBorn}
\end{figure}

In all processes we evaluate the Feynman diagrams, performing the Dirac traces with the TRACER package \cite{Jamin:1991dp}, and we integrate over the virtual loop and real emission momenta. Ghost graphs are taken into account to simplify the polarization sums of external gluons. Real contributions are computed using traditional partial fractioning and master integrals, while virtual loop contributions are calculated using the FeynCalc \cite{Shtabovenko:2016sxi} and Feynhelpers \cite{Shtabovenko:2016whf} packages. Dimensional regularization is used to regulate infrared and ultraviolet singularities, working in a space-time of $d=4-2 \epsilon $ dimensions. As it will be addressed in the following subsections, dimensional regularization requires a special treatment for the Levi-Civita tensor $\epsilon^{\mu \nu \rho \sigma}$ and the Dirac matrix $\gamma^5$, which is present in both the electroweak couplings and quirality projectors. The factorization of the mass singularities in the polarized case is also affected by the scheme used to deal with $\gamma^5$.

\subsubsection{Treatment of the matrix $\gamma_5$ and axial couplings}

While dimensional regularization involves working in $d$-dimensional space-time, the $\gamma^5$ matrix and the $\epsilon^{\mu \nu \rho \sigma}$ tensor are only well defined in the four-dimensional space-time. A consistent way to treat $\gamma^5$ and $\epsilon^{\mu \nu \rho \sigma}$ in $d$-dimensions is the HVBM scheme \cite{tHooft:1972tcz,Breitenlohner:1977hr}, which splits the $d$-dimensional Minkowsky space into the usual four-dimensional one and a $(d-4)$-dimensional subspace where, for instance, the $(d-4)$-dimensional part of the $\gamma^{\mu}$ matrices, represented as $\hat{\gamma}^{\mu}$, commutes with the strictly four-dimensional $\gamma^5$. Calculations in the HVBM scheme are algebraically more involved because for unobserved momenta $p$ and $\gamma^{\mu}$ matrices one needs to take into account both their four-dimensional part ($\tilde{p}$, $\tilde{\gamma^{\mu}}$) and their $(d-4)$-dimensional one ($\hat{p}$, $\hat{\gamma^{\mu}}$), which follow different algebraic treatment. Integration over `hat' momenta $\hat{p}$ in the real contributions is performed as established in  Ref.\cite{Gordon:1993qc}.

One of the first issues related to the HVBM scheme is the definition of the electroweak vertices. The standard definition of the boson vertices, taking into account the flavours $f_1$ and $f_2$ of the involved quarks, can be expressed as

\begin{equation}
    -i e \gamma_{\mu} \left( L_{f_2 f_1} \frac{1-\gamma_5}{2} + R_{f_2 f_1} \frac{1+\gamma_5}{2} \right),
\label{eq_vert4D}
\end{equation}

\noindent where the left and right-handed coupling are given by the following expressions, depending on the type of boson

\begin{equation}
    \begin{aligned}[c]
        W^{-}: & \quad L_{f_2 f_1}=\frac{1}{\sqrt{2} \sin \theta_W}(\tau_{+})_{f_2 f_1} U_{f_2 f_1}, \quad R_{f_2 f_1}=0, \\
        W^{+}: & \quad L_{f_2 f_1}=\frac{1}{\sqrt{2} \sin \theta_W}(\tau_{-})_{f_2 f_1} U^{\dagger}_{f_2 f_1}, \quad R_{f_2 f_1}=0 \\
        Z^{0}: & \quad L_{f_2 f_1}=\frac{1}{\sin (2 \theta_W)}(\tau_3)_{f_2 f_1}- \delta_{f_2 f_1} e_{f_1} \tan \theta_W, \quad R_{f_2 f_1}=-\delta_{f_2 f_1} e_{f_1} \tan \theta_W, \\
        \gamma: & \quad L_{f_2 f_1}=\delta_{f_2 f_1} e_{f_1}, \quad R_{f_2 f_1}=\delta_{f_2 f_1} e_{f_1}. \\
    \end{aligned}
\end{equation}

\noindent Here $\theta_W$ is the electroweak mixing angle, $\tau_{\pm}=(\tau_1 \pm i \tau_2)/2$ and $\tau_3$ are the weak isospin Pauli matrices, $U$ is the CKM mixing matrix, and $e_f$ is the electric charge of the corresponding quark ($e_f=\frac{2}{3}$ for $u,c,t$ and $e_f=-\frac{1}{3}$ for $d,s,b$). However, since the commutation properties of the combination $\gamma_{\mu} \gamma_{5}$ depends on whether $\mu \leq 3$ or $\mu > 3$, a symmetric definition of the vertices is needed, which corresponds to the replacement $\gamma_{\mu}(1-\gamma_5) \rightarrow \frac{1}{2}(1+\gamma_5)\gamma_{\mu}(1-\gamma_5)$ \cite{Korner:1989is,Buras:1989xd}. After this substitution, the vertex expression of Eq.(\ref{eq_vert4D}) in $d$-dimensions can be rewritten as

\begin{equation}
    -i e \tilde{\gamma}_{\mu} \left( L_{f_2 f_1} \frac{1-\gamma_5}{2} + R_{f_2 f_1} \frac{1+\gamma_5}{2} \right),
\end{equation}

\noindent where the $(d-4)$-dimensional part of $\gamma_{\mu}$ cancels within the left and right-handed terms.

In the HVBM scheme the presence of $\gamma_5$ in both the symmetrized vertex and the quirality projectors over initial partons lead to anomalous terms of order $\mathcal{O}(d-4)$, which in combination with divergent terms can result in spurious finite contributions. Infrared anomalous terms cancel out between virtual and real contributions \cite{Korner:1985uj}, but in the case of the renormalization of ultraviolet divergencies and mass factorization appropriate additional counterterms are required.

\subsubsection{Renormalization and the axial vertex anomaly}

For the renormalization of the virtual contributions to the NLO we adopt de usual $\overline{MS}$ scheme. In this scheme the running strong coupling constant $\alpha_s(\mu^2)$ at a scale $\mu^2$ satisfies the corresponding renormalization-group equation

\begin{equation}
    \mu^2 \frac{d}{d \mu^2} \alpha_s (\mu^2) =-\alpha_s (\mu^2) \left[ \beta_0 \frac{\alpha_s (\mu^2)}{2 \pi} + \beta_1 \left( \frac{\alpha_s (\mu^2)}{2 \pi} \right)^2+ \mathcal{O}(\alpha_s^3) \right],
\end{equation}

\noindent with the $\beta_i$ factors given by

\begin{equation}
    \begin{aligned}[c]
        \beta_0 &= \frac{11}{6} C_A-\frac{1}{3} \sum_f (1), \\
	\beta_1 &= \frac{17}{6} C_A^2-\left(\frac{5}{6} C_A + \frac{1}{2} C_F \right) \sum_f (1),
    \end{aligned}
\end{equation}

\noindent where the sum $\sum_f$ runs over all the flavours below the mass threshold.

However, within the HVBM scheme, since $\gamma^5$ no longer anticommutes with $\gamma^{\mu}$ an additional finite renormalization of the axial quark current needs to be taken into account at NLO \cite{Larin:1993tq}. The renormalization constant related to the axial current $Z_5$ is no longer equal to one. The corresponding one loop expression is

\begin{equation}
    Z_5 =1-\frac{\alpha_s}{2 \pi} \ 2 C_F+\mathcal{O}(\alpha_s^2).
\end{equation}

\noindent This term gives rise to the following additional finite counterterms that have to be added to the $q\overline{q}$ and $qg$ channels in order cancel the spurious axial terms appearing in the virtual contributions
\begin{equation}
    \begin{aligned}[c]
      (\Delta)  C^{q\overline{q}}=&-\frac{\alpha_s}{2 \pi} \left( 2 \ C_F\right) \frac{d (\Delta) \sigma^{axial}_{q\overline{q}}}{d^3 Q}, \\
	(\Delta)C^{qg}=&-\frac{\alpha_s}{2 \pi} \left( 2 \ C_F\right) \frac{d (\Delta) \sigma^{axial}_{qg}}{d^3 Q},
    \end{aligned}
\label{eq_axrenct}
\end{equation}

\noindent where $d (\Delta) \sigma^{axial}_{ij}/d^3 Q$ is the axial part of the (polarized) unpolarized Born cross section corresponding to either the $q \overline{q} \rightarrow V g$ process or the  $q g \rightarrow V q$ Compton process. This axial cross section is obtained by considering only the terms depending on the axial coupling to the boson $V$.

\subsubsection{Factorization}

Mass singularities appear due to the collinear emission of massless partons from one of the incoming partons. These singular terms are detached at a scale $\mu_F^2$ to be incorporated into the (polarized) unpolarized distribution functions $(\Delta)f_a^h(x,\mu_F^2)$. For the unpolarized case, the results are calculated in the usual $\overline{MS}$ scheme. However, in the polarized case we use the conventional variation of the $\overline{MS}$ scheme which takes into account some helicity-conservation violations that arise in the HVBM scheme, particularly in the quark-antiquark annihilation process.

\begin{figure}
\centering
  \includegraphics[width=0.4\linewidth]{Fig_Factorization.pdf}
  \caption{Schematic of the parton emission processes.}
  \label{fig_fact}
\end{figure}

The factorization subtractions \footnote{Subtractions that need to be applied to the corresponding partonic channel.} due to collinear emission from an initial parton $a$, as illustrated in Fig.\ref{fig_fact} for gluon emission, are given by the following convolution

\begin{equation}
    d (\Delta) \hat{\sigma}^{fact}_{a b \rightarrow V X}=\frac{\alpha_s}{2 \pi} \left\{ \left[ -\frac{1}{\hat{\epsilon}} (\Delta) P_{ca}(z) \left(\frac{\mu^2}{\mu_F^2} \right)^{\epsilon}+ (\Delta) f_{ca} (z)\right] \circledast d (\Delta) \sigma_{cb \rightarrow V X} \right\},
\end{equation}

\noindent where $\mu$ is the scale of dimensional regularization and $\mu_F$ represents the factorization scale. Here $1/\hat{\epsilon}=1/\epsilon-\gamma_E+\ln 4 \pi$ as in the $\overline{MS}$ scheme, $(\Delta) P_{ca}(z)$ is the (polarized) unpolarized Altarelli-Parisi splitting function of a parton $a$ emitting a parton $c$ carrying momentum fraction $z$, $d (\Delta) \sigma_{cb \rightarrow V X}$ is the (polarized) unpolarized $d$-dimensional Born cross section of the partonic process $c b \rightarrow V X$, and $(\Delta) f_{ab}(z)$ represents the freedom on subtracting additional finite terms. In the unpolarized case we rely on the usual $\overline{MS}$ scheme, which corresponds to $f_{ij}(z)=0$, but in the polarized case when the subtraction involves emission of quarks by other quarks, when working with the $\Delta P_{qq}(z)$ kernel, we use $\Delta f_{qq}(z)=-4 \ C_F (1-z)$. This is the particular scheme is used in the definition of the polarized parton distribution sets \cite{deFlorian:2009vb}.

\section{NLO corrections}
\label{sect2}

In this section we present the NLO contributions to the polarized cross section of electroweak boson production at large transverse momentum in hadronic collisions. For the sake of completeness we also provide the unpolarized results, which are in agreement with those of Gonsalves et al. \cite{Gonsalves:1989ar}, with the distinction that we corroborated the calculation fully within the consistent HVBM scheme for the axial coupling. The polarized cross section for the virtual photon in the non-singlet case partially agrees with the previous results of Ref.\cite{Chang:1997ik} due to an overlooked sign in one of the interferences.

\subsection{Quark-antiquark annihilation and scattering}

The formulas of the inclusive cross section are presented in the same way as in Ref.\cite{Gonsalves:1989ar}, with the expression separated in contributions with different axial-vector and/or flavour structure. The quark-antiquark cross section formula is given by

\begingroup
\small
\begin{flalign}
 E_Q \frac{d (\Delta) \sigma^{q\overline{q}}}{d^3 Q}=&\frac{\alpha \ \alpha_s^{\overline{MS}}(\mu^2) \ C_F}{s \ N_c} \Bigg\{ \delta (s_{23}) \ (\Delta) B_{q\overline{q}}(s,t,u,Q^2) \left( |L_{21}|^2+|R_{21}|^2 \right)+\frac{\alpha_s^{\overline{MS}}(\mu^2)}{2 \pi} \notag & \\
 & \times \Big\{ \Big[ \delta (s_{23}) \big[ (\Delta) V^{(3)}_{q\overline{q}}(s,t,u,Q^2) + (\Delta) V^{(1)}_{q\overline{q}}(s,t,u,Q^2) \sum_f (1) \big] \notag \\
 &\quad + (\Delta) G_{q\overline{q}}(s,t,u,Q^2) + (\Delta) F_{aa}(s,t,u,Q^2) \sum_f (1) \Big] \left( |L_{21}|^2+|R_{21}|^2 \right) \notag \\
 &\quad + \Big[ \delta (s_{23}) \ (\Delta) V^{(2)}_{q\overline{q}}(s,t,u,Q^2) + (\Delta) F_{ab}(s,t,u,Q^2) \Big] \delta_{12} \left(L_{11}-R_{11}\right) \sum_f \left(L_{ff}-R_{ff}\right) \notag \\
 &\quad + (\Delta) F_{bb}(s,t,u,Q^2) \delta_{12} \sum_f \sum_{f'} \left( |L_{ff'}|^2+|R_{ff'}|^2 \right) \notag \\
 &\quad +\Big[ (\Delta) F_{ac}(s,t,u,Q^2) + (\Delta) F_{ad}(s,t,u,Q^2) \Big] \left( |L_{21}|^2+|R_{21}|^2 \right) \notag \\
 &\quad +\Big[ \delta_{12} (\Delta) F_{bc}(s,t,u,Q^2) + (\Delta) F_{cc}(s,t,u,Q^2) \Big] \sum_f \left( |L_{f1}|^2+|R_{f1}|^2 \right) \notag \\
 &\quad +\Big[ \delta_{12} (\Delta) F_{bd}(s,t,u,Q^2) + (\Delta) F_{dd}(s,t,u,Q^2) \Big] \sum_f \left( |L_{2f}|^2 +|R_{2f}|^2 \right) \notag \\
 &\quad + (\Delta) F^{LL}_{cd}(s,t,u,Q^2) \left( L_{11} L_{22}+R_{11} R_{22} \right) \notag \\
 &\quad + (\Delta) F^{LR}_{cd}(s,t,u,Q^2) \left( L_{11} R_{22}+R_{11} L_{22} \right)\Big\} \Bigg\}.
\label{eq_masterqqb}
\end{flalign}
\endgroup

\noindent Here the $(\Delta) B_{q\overline{q}}(s,t,u,Q^2)$ function corresponds to the (polarized) unpolarized Born diagrams contributions \footnote[1]{The axial Born cross section needed for the axial renormalization counterterm in Eq.(\ref{eq_axrenct}) is given by the expression $E_Q \frac{d (\Delta) \sigma^{axial}_{q\overline{q}}}{d^3 Q}=\frac{\alpha \ \alpha_s^{\overline{MS}}(\mu^2) \ C_F}{s \ N_c} (\Delta) B_{q\overline{q}}(s,t,u,Q^2) \frac{\left( L_{21}-R_{21} \right)^2}{2}$. The expression is analogous in the $qg$ case.}

\small
\begin{equation}
 \Delta B_{q\overline{q}}(s,t,u,Q^2)=-\frac{t^2+u^2+2 s Q^2}{t u}, \quad B_{q\overline{q}}(s,t,u,Q^2)=-\Delta B_{q\overline{q}}(s,t,u,Q^2).
\label{eq_Bornqqb}
\end{equation}
\normalsize

\noindent The rest of the functions are given in the Appendix. The $(\Delta) V^{(i)}_{q\overline{q}}(s,t,u,Q^2)$ functions group all the NLO contributions proportional to $\delta(s_{23})$ arising mainly from virtual diagrams of Fig.\ref{fig_FeynVirt}, but also from infrared emission from the real ones. The $(\Delta) G_{q\overline{q}}(s,t,u,Q^2)$ are the contributions from diagrams with two final gluons of Fig.\ref{fig_FeynqqGG}, while $(\Delta) F(s,t,u,Q^2)$ functions represent the contributions from the diagrams of Fig.\ref{fig_Feynqqbqqb}, with subscripts indicating the particular groups of diagrams involved in each subprocess. Most of these functions differ by a sign between the polarized and unpolarized case, with the only exceptions being $(\Delta) F_{cc}(s,t,u,Q^2)$, $(\Delta) F_{dd}(s,t,u,Q^2)$ and $(\Delta) F^{LR}_{cd}(s,t,u,Q^2)$. This is related to helicity conservation for massless quarks. The $F^{LR}_{cd}(s,t,u,Q^2)$ function is actually equal to its polarized counterpart.

\begin{figure}
\centering
  \includegraphics[width=0.6\linewidth]{Fig_Feynman_Virtual.pdf}
  \caption{Diagrams corresponding to the virtual corrections to the Born $q \overline{q} \rightarrow V g$ process. The $V$ functions in the cross section formula correspond to these diagrams.}
  \label{fig_FeynVirt}
\end{figure}

\begin{figure}
\centering
  \centering
  \includegraphics[width=0.6\linewidth]{Fig_Feynman_qq-gg.pdf}
  \caption{Diagrams which contribute to the annihilation process $q \overline{q} \rightarrow V g g$. Ghost graphs are not depicted. The $G$ functions in the cross section formula correspond to these diagrams.}
  \label{fig_FeynqqGG}
\end{figure}

\begin{figure}
  \centering
  \includegraphics[width=0.6\linewidth]{Fig_Feynman_qqb-qqb.pdf}
  \caption{Diagrams which contribute to the process $q \overline{q} \rightarrow V q \overline{q}$. The $F$ functions in the cross section formula correspond to these diagrams.}
  \label{fig_Feynqqbqqb}
\end{figure}

The cross section of the process $\overline{q} q \rightarrow V X$ can be obtained from Eq.(\ref{eq_masterqqb}) by substituting the corresponding couplings
\begingroup
\small
\begin{flalign}
    \overline{q} q \rightarrow V X: \quad L \leftrightarrow -R^{\dagger}. &&
\end{flalign}
\endgroup

\subsection{Quark-gluon Compton scattering}

The cross section of the Compton process is given by
\begingroup
\small
\begin{flalign}
 E_Q \frac{d (\Delta) \sigma^{qg}}{d^3 Q}=\frac{\alpha \ \alpha_s^{\overline{MS}}(\mu^2) \ C_F}{s \ (N_c^2-1)} & \Bigg\{ \delta (s_{23}) \ (\Delta) B_{qg}(s,t,u,Q^2) \sum_f \left( |L_{f1}|^2+|R_{f1}|^2 \right) \notag & \\
 +& \frac{\alpha_s^{\overline{MS}}(\mu^2)}{2 \pi} \Big\{ \Big[ \delta (s_{23}) \big[ (\Delta) V^{(3)}_{qg}(s,t,u,Q^2) + (\Delta) V^{(1)}_{qg}(s,t,u,Q^2) \sum_f (1) \big] \notag \\
 &\quad \quad + (\Delta) G_{qg}(s,t,u,Q^2) \Big] \sum_f \left( |L_{f1}|^2+|R_{f1}|^2 \right) \notag \\
 & \quad + \delta (s_{23}) (\Delta) V^{(2)}_{qg}(s,t,u,Q^2) \left(L_{11}-R_{11}\right) \sum_f \left(L_{ff}-R_{ff}\right) \Big\} \Bigg\},
\label{eq_masterqg}
\end{flalign}
\endgroup
\noindent where in this case the (polarized) unpolarized Born terms function reads
\small
\begin{equation}
 \Delta B_{qg}(s,t,u,Q^2)=\frac{(s-t) \left(Q^2+u\right)}{s t}, \quad B_{qg}(s,t,u,Q^2)=-\frac{s^2+t^2+2 u Q^2}{s t}.
\label{eq_Bornqg}
\end{equation}
\normalsize
\noindent The $(\Delta) G_{qg}(s,t,u,Q^2)$ terms come from diagrams obtained by crossing from those in Fig.\ref{fig_FeynqqGG}, while the $(\Delta) V^{(i)}_{qg}(s,t,u,Q^2)$ are again the $\delta(s_{23})$ contributions originating from virtual diagrams contributions, which can be obtained from Fig.\ref{fig_FeynVirt} by crossing, and from the same real contributions. Their explicit formulas are given in the Appendix.

The $qg$ process presents the greatest difference between the polarized and unpolarized cross sections. As a matter of fact, Eq.(\ref{eq_Bornqg}) reveals that at the Born level they have two distinct behaviours: one is symmetric in the variables $s$ and $t$, while the other one is completely antisymmetric in them.

The three variations of the Compton process can be obtained from Eq.(\ref{eq_masterqg}) with the following substitutions
\begingroup
\small
\begin{flalign}
    \overline{q} g \rightarrow V X:& \quad L \leftrightarrow -R^{\dagger}, & \notag \\
    g q \rightarrow V X:& \quad t \leftrightarrow u, \quad f_1 \leftrightarrow f_2, & \notag \\
    \overline{q} g \rightarrow V X:& \quad t \leftrightarrow u, \quad f_1 \leftrightarrow f_2, \quad L \leftrightarrow -R^{\dagger}.&
\end{flalign}
\endgroup

\subsection{Gluon-gluon fusion}

The cross section of the gluon-gluon fusion process, which only contributes at NLO, is simply given by

\begingroup
\small
\begin{flalign}
 E_Q \frac{d (\Delta) \sigma^{gg}}{d^3 Q}&=\frac{\alpha \ \alpha_s^{\overline{MS}}(\mu^2) \ N_c \ C_F}{s \ (N_c^2-1)^2} \left( \frac{\alpha_s^{\overline{MS}}(\mu^2)}{2 \pi} \right) (\Delta) G_{gg}(s,t,u,Q^2) \sum_f \sum_{f'} \left( |L_{ff'}|^2+|R_{ff'}|^2 \right), &
\end{flalign}
\label{eq_mastergg}
\endgroup

\noindent where the $(\Delta) G_{gg}(s,t,u,Q^2)$ functions, corresponding to the diagrams obtained by crossing from the ones in Fig.\ref{fig_FeynqqGG}, are given in the Appendix. The $G_{gg}(s,t,u,Q^2)$ and $\Delta G_{gg}(s,t,u,Q^2)$ functions are fairly similar apart from a sign and a few additional terms.

\subsection{Quark-quark scattering}

Finally, the last NLO contribution is the quark-quark initiated, whose cross section is given by

\begingroup
\small
\begin{flalign}
 E_Q \frac{d (\Delta) \sigma^{qq}}{d^3 Q}=&\frac{\alpha \ \alpha_s^{\overline{MS}}(\mu^2) \ C_F}{s \ N_c} \left( \frac{\alpha_s^{\overline{MS}}(\mu^2)}{2 \pi} \right) \left( \frac{1}{2} \right) \notag & \\
 & \times \Big\{ \big[ (\Delta) H_{aa}(s,t,u,Q^2) + (\Delta) H_{cc}(s,t,u,Q^2) \big] \sum_f \left( |L_{f1}|^2+|R_{f1}|^2 \right) \notag \\
 &+ \big[ (\Delta) H_{bb}(s,t,u,Q^2) + (\Delta) H_{dd}(s,t,u,Q^2) \big] \sum_f \left( |L_{2f}|^2+|R_{2f}|^2 \right) \notag \\
 &+ (\Delta) H_{ac}(s,t,u,Q^2) \left( |L_{21}|^2+|R_{21}|^2 \right) + (\Delta) H_{bd}(s,t,u,Q^2) \left( |L_{12}|^2+|R_{12}|^2 \right) \notag \\
 &+\Big[ (\Delta) H_{ad}(s,t,u,Q^2) + (\Delta) H_{bc}(s,t,u,Q^2) \Big] \delta_{12} \sum_f \left( |L_{f1}|^2+|R_{f1}|^2 \right) \notag \\
 &+ \Big[ (\Delta) H^{LL}_{ab}(s,t,u,Q^2) + (\Delta) H^{LL}_{cd}(s,t,u,Q^2) \Big] \left( L_{11} L_{22}+R_{11} R_{22} \right) \notag \\
 &+ \Big[ (\Delta) H^{LR}_{ab}(s,t,u,Q^2) + (\Delta) H^{LR}_{cd}(s,t,u,Q^2) \Big] \left( L_{11} R_{22}+R_{11} L_{22} \right) \Big\},
\label{eq_masterqq}
\end{flalign}
\endgroup

\noindent where the $(\Delta) H(s,t,u,Q^2)$ functions correspond to the contributions from the diagrams in Fig.\ref{fig_Feynqq}. The global $1/2$ factor is a statistical factor needed to avoid double counting events due to the full phase space integration and quantum number summation over the two quarks in the final state. Most of the $H$ functions are related to the $F$ ones from the quark-antiquark process. The $H_{ac}(s,t,u,Q^2)$, $H_{bd}(s,t,u,Q^2)$, $H_{ad}(s,t,u,Q^2)$, $H_{bc}(s,t,u,Q^2)$, $H^{LR}_{ab}(s,t,u,Q^2)$ and $H^{LR}_{cd}(s,t,u,Q^2)$ functions are equal to their polarized counterparts, while the $H^{LL}_{ab}(s,t,u,Q^2)$ and $H^{LL}_{cd}(s,t,u,Q^2)$ differ only by a sign. The $H_{bb}(s,t,u,Q^2)$ and $H_{dd}(s,t,u,Q^2)$ functions, as it is the case with $F_{cc}(s,t,u,Q^2)$ and $F_{dd}(s,t,u,Q^2)$, have significant differences with their polarized counterparts.

\begin{figure}
\centering
  \includegraphics[width=0.6\linewidth]{Fig_Feynman_qq-qq.pdf}
  \caption{Diagrams corresponding to the $q q \rightarrow V q q$ process. The $H$ functions in the cross section formula correspond to these diagrams.}
  \label{fig_Feynqq}
\end{figure}

The cross section of antiquark-antiquark annihilation is obtained from Eq.(\ref{eq_masterqq}) by the substitution
\begingroup
\small
\begin{flalign}
    \overline{q} \overline{q} \rightarrow V X: \quad L \leftrightarrow -R^{\dagger}. &&
\end{flalign}
\endgroup

With the $q \overline{q}$ and $q q$ cross sections it is possible to calculate the non-singlet cross section, which is given by
\small
\begin{equation}
 E_Q \frac{d (\Delta) \sigma^{NS}}{d^3 Q}=E_Q \frac{d (\Delta) \sigma^{q \overline{q}}}{d^3 Q}-E_Q \frac{d (\Delta) \sigma^{q q}}{d^3 Q}.
\end{equation}
\normalsize

The non-singlet polarized cross section for virtual photon production has been calculated by Field et al. in Ref.\cite{Chang:1997ik}. However, our results do not fully agree with theirs because they assumed that the $\Delta H^{LL}_{ab}(s,t,u,Q^2)$, $\Delta H^{LL}_{cd}(s,t,u,Q^2)$, $\Delta H^{LR}_{ab}(s,t,u,Q^2)$ and $\Delta H^{LR}_{cd}(s,t,u,Q^2)$ contributions of the $q q$ process cancel against the corresponding $\Delta F^{LL}_{cd}(s,t,u,Q^2)$ and $\Delta F^{LR}_{cd}(s,t,u,Q^2)$ functions of the $q \overline{q}$ process, which is not correct since they actually differ by a sign and add up to an additional contribution. This error was carried over from a mistake in the original unpolarized EMP calculation \cite{Ellis:1981hk}, which it is pointed out in Ref.\cite{Gonsalves:1989ar}.

\section{Phenomenological Results}
\label{sect3}

In this section we compute the corresponding cross sections to evaluate the impact of the polarized NLO corrections. We apply them to  RHIC kinematics, that is, polarized proton-proton collisions at a center-of-mass energy of $\sqrt{S}=510$ GeV. We use the MMHT2014nlo \cite{Harland-Lang:2014zoa} and DSSV \cite{deFlorian:2009vb} PDFs sets for the unpolarized and polarized processes, respectively. For the weak gauge boson we use the masses $M_W=80.379$ GeV and $M_Z=91.1876$ GeV, with an electromagnetic coupling constant $\alpha=1/127.918$ and the Weinberg angle given by $\sin{\theta_W}^2=0.2312$. The CKM matrix is constructed from the values $|V_{ud}|=0.9742$, $|V_{us}|=0.2243$, $|V_{cd}|=0.218$ and $|V_{cs}|=0.996$ using the universality property. The number of active quark flavors are limited by the thresholds set by their mass, namely $m_c=1.275$ GeV and $m_b=4.18$ GeV. The QCD coupling $\alpha_s$ is evaluated at NLO with the same quark thresholds, with a value of $\alpha_s(M_Z)=0.127$.

\begin{figure}[h]
\centering
  \includegraphics[width=0.5\linewidth]{Fig_plotZuncert.pdf}
  \caption{On-shell $Z$ production cross section in $pp$ collisions at $510$ GeV. Both LO and NLO are presented with their respective confidence bands.}
  \label{fig_plotZunc}
\end{figure}

We begin by analysing the reduction of theoretical uncertainties. These are estimated by studying the dependence of the complete NLO results on the renormalization and factorization scales due to the truncation of the perturbative series expansion. In Fig.\ref{fig_plotZunc} we present the polarized $q_T$ distributions of on-shell $Z$ boson production at both LO and NLO accuracy at $y=0$, with their respective confidence bands. Since in this process we have two physical scales given by $Q$ and $q_T$, it is convenient to define a physical scale in between them given by $Q_s^2=(Q^2+q_T^2)/4$. The central curves are obtained by fixing the renormalization and factorization scales to $\mu_R=\mu_F=Q_s$. The bands are obtained by varying $\mu_R$ and $\mu_F$ simultaneously and independently in the range $0.5 Q_s \ \leq \ \mu_R,\mu_F \ \leq \ 2 Q_s$ with the constraint $0.5 \leq \mu_F/\mu_R \leq 2$. As expected in perturbative QCD, the second order corrections yield a significant contribution to the Born cross section while reducing uncertainties considerabily.

For a more detailed analysis on the impact of the two scales, in Fig.\ref{fig_plotZscldep} we study the scale dependence of the cross section at a fixed value of $q_T$. The LO curves show stronger dependence on the scale. In this case the $\mu_R$ dependence corresponds simply to the variations of $\alpha_s$ with the scale, while $\mu_F$ dependence arises exclusively from the PDFs. At NLO these dependences are corrected by logarithmic terms proportional to the $\beta$ function and the splitting functions, respectively, yielding more stable results. We do not show the results for other particular values of $q_T$ since they are fairly similar. Throughout the rest of this section we will set the two scales equal to $Q_s$.

\begin{figure}[!b]
\centering
  \includegraphics[width=0.5\linewidth]{Fig_plotZscaledep.pdf}
  \caption{Cross section dependence on the renormalization and factorization scales. The green curves are for $\mu=\mu_R$ and $\mu_F=Q_s$, the red ones for $\mu=\mu_F$ and $\mu_R=Q_s$, and the black ones for $\mu=\mu_F=\mu_R$.}
  \label{fig_plotZscldep}
\end{figure}

The size of the higher order QCD corrections to hadronic processes are usually presented in terms of the `$K$-factor', the ratio between the NLO results over the LO ones. However, in order to make this ratio meaningful in the polarized case, the same NLO-evolved parton densities are used to calculate both NLO and LO quantities \cite{deFlorian:1998qp}. This is due to the low constraints to polarized PDFs available from experimental data, especially to $\Delta g$, which can in some cases give very different results when fits are performed at LO or at NLO. For this reason, polarized $K$-factors can get artificially large or small when the gluons are involved.

In Fig.\ref{fig_plotkf} we present the $K$-factors for the production of on-shell weak bosons at fixed rapidity $y=0$ and $y=1$, for both the polarized and unpolarized cases. The polarized $K$ factors reveal a larger contribution of the NLO corrections for $W$ production, and they show a stronger variation along $q_T$ compared to their unpolarized counterparts. The polarized $Z$ production present a remarkably lower $K$ factor at larger $q_T$. At rapidity $y=1$ the $K$-factors are slightly higher, specially at low $q_T$. In this case the bosons $W^{+}$ and $W^{-}$ show similar behaviours.

\begin{figure}
\begin{subfigure}{0.49\textwidth}
\includegraphics[width=\linewidth]{Fig_plotkfy=0.pdf}
\end{subfigure}
\hspace*{\fill}
\begin{subfigure}{0.49\textwidth}
\includegraphics[width=\linewidth]{Fig_plotkfy=1.pdf}
\end{subfigure}
\caption{NLO $K$ factors for on-shell weak boson $V$ production with rapidities $y=0$ and $y=1$. The solid and dashed curves correspond to the polarized and unpolarized processes, respectively.} \label{fig_plotkf}
\end{figure}

The $q_T$ dependence of the polarized $K$ factors is related to the different behaviour of the particular sub-channels contributing to the process. For instance, the changes in the concavity in the curves for $y=0$ (barely visible in the $Z$ boson case) are related to the sign change  of the $qg$-channel contribution. This can be observed in Fig.\ref{fig_plotchann} for $W^{+}$ production, where we present the sub-channel polarized cross section ratios. In all cases the cross section is governed by the quark initiated channels at low $q_T$, but the $qg$-channel becomes the most relevant one at high $q_T$. The $q_T$ range at which this dominance begins depends on the particular boson studied and its rapidity, taking place at lower $q_T$ values in the case $y=1$. The $gg$ contribution, that only arises at order $\mathcal{O}(\alpha_s^2)$, is always negligible at the energy of this process.

\begin{figure}[h]
\begin{subfigure}{0.49\textwidth}
\includegraphics[width=\linewidth]{Fig_plotcanalesy0.pdf}
\end{subfigure}
\hspace*{\fill}
\begin{subfigure}{0.49\textwidth}
\includegraphics[width=\linewidth]{Fig_plotcanalesy1.pdf}
\end{subfigure}
\caption{Sub-channel contributions to on-shell $W^{+}$ boson production with rapidities $y=0$ and $y=1$. Both NLO and LO decompositions are given.} \label{fig_plotchann}
\end{figure}

Finally in Fig.\ref{fig_plotall} we show the longitudinal double spin asymmetries for weak boson production at different rapidities, which is defined as

\begin{equation}
 A_{LL}=\frac{d \Delta \sigma}{d \sigma}.
\end{equation}

\noindent The asymmetries grow with $q_T$, reaching sizeable values, and are larger at higher rapidity. The NLO corrections to the asymmetries are generally low, especially at rapidity $y=1$. A similar behaviour has been observed in many other processes, where QCD corrections tend to compensate in the ratios of polarized and unpolarized cross sections. However, in the case of $Z$ production the NLO corrections in the polarized case are clearly bigger than the unpolarized ones, accounting for lower values of asymmetry at large $q_T$. This can already be seen in Fig.\ref{fig_plotkf} since the ratios of the asymmetries at the two orders are equal to the ratio of the respective $K$ factors, following the relation $A^{NLO}=A^{LO} \ \Delta K/K$, and the K factors of $Z$ production are the ones that show larger discrepancies between the polarized and unpolarized cases, especially at $y=0$.

\begin{figure}[h]
\begin{subfigure}{0.49\textwidth}
\includegraphics[width=\linewidth]{Fig_plotally=0.pdf}
\end{subfigure}
\hspace*{\fill}
\begin{subfigure}{0.49\textwidth}
\includegraphics[width=\linewidth]{Fig_plotally=1.pdf}
\end{subfigure}
\caption{Longitudinal double spin asymmetries in on-shell weak boson $V$ production at rapidities $y=0$ and $y=1$.} \label{fig_plotall}
\end{figure}

\section{Summary}
\label{sec:summa}

In this work, we have presented the first complete calculation at next-to-leading order in perturbative QCD of the cross section for the production of electroweak gauge bosons with large transverse momentum $q_T$ in double polarized hadronic collisions. The calculation was done fully in the consistent HVBM scheme for the treatment of $\gamma^5$ and the tensor $\epsilon_{\mu \nu \rho \sigma}$ in dimensional regularization, carefully dealing with issues arising from the $(d-4)$-dimensional part of the momenta and $\gamma^{\mu}$ matrices. We also kept the unpolarized results, easily obtained as a by-product of the polarized calculation, corroborating the previous well-known results \cite{Gonsalves:1989ar}, but this time fully within the HVBM scheme. The polarized results are also in agreement with a previous partial calculation of the non-singlet polarized production of virtual photons \cite{Chang:1997ik}, aside from a missed contribution in their work.

Using our results, we studied in some detail the phenomenological impact of the polarized NLO corrections in weak boson production at RHIC, that is, polarized $pp$ collisions with a center-of-mass energy of 510 GeV. We analysed the scale dependence in polarized $Z$ boson production to assess the reliability of the theoretical predictions. We found that the NLO correction have a greatly reduced scale dependence with respect to the Born results. We also studied the polarized $K$ factors, the ratio between the NLO and LO results, comparing them to their unpolarized counterparts. The polarized $K$ factors show larger variations and a more complex dependence in $q_T$, having a slightly higher value than the unpolarized ones except for the $Z$ boson case. The effect of NLO corrections on the double spin asymmetries is generally low for $W$ production, but visibly affects the $Z$ production at higher $q_T$.

\section*{Acknowledgements}
\label{sec:akno}

We are grateful to Werner Vogelsang for helpful discussions. This work has been partially supported by Conicet and ANPCyT.