# Explaining Reinforcement Learning Policies through Counterfactual
Trajectories

Julius Frost  Olivia Watkins  Eric Weiner  Pieter Abbeel  Trevor Darrell
Bryan Plummer  Kate Saenko

###### Abstract

In order for humans to confidently decide where to employ RL agents for real-
world tasks, a human developer must validate that the agent will perform well
at test-time. Some policy interpretability methods facilitate this by
capturing the policy’s decision making in a set of agent rollouts. However,
even the most informative trajectories of training time behavior may give
little insight into the agent’s behavior out of distribution. In contrast, our
method conveys how the agent performs under distribution shifts by showing the
agent’s behavior across a wider trajectory distribution. We generate these
trajectories by guiding the agent to more diverse unseen states and showing
the agent’s behavior there. In a user study, we demonstrate that our method
enables users to score better than baseline methods on one of two agent
validation tasks.

Explainable AI, Reinforcement Learning

##  1 Introduction

In most of the benchmarks where reinforcement learning (RL) has shown
impressive results (Bellemare et al., 2013), (Tassa et al., 2018), agents are
trained and tested in the same environment. In contrast, many practical
applications involve sim-to-real transfer or a real-world train/test
distribution shift which can hurt test-time performance (Eysenbach et al.,
2020), (Akkaya et al., 2019), (Rusu et al., 2017), (Peng et al., 2018).
Therefore it important to validate the performance and behavior of RL agents
in the presence of a distribution shift. To do this, we take a step beyond
summary metrics and focus on explainability techniques which help humans gain
a more in-depth understanding of agent behavior.

Figure 1: Our method seeks out more diverse states from which to show rollouts
of agent behavior to the user.

Previous work in explainability often aims to help users understand and
predict model behavior by explaining a set of model predictions (Hoffman et
al., 2018). In the reinforcement learning setting, it may be desirable to
select the most informative segments that show the behavior of the agent. With
this motivation, critical states methods show the agent’s behavior in states
which the agent thinks are important (Huang et al., 2018), (Amir & Amir,
2018). However, even if these methods successfully explain policy behavior, we
cannot conclude if the behavior generalizes when the distribution shifts in
the environment.

Our work builds upon the core idea of creating a set of informative
trajectories of the agent’s behavior, but we focus on showing trajectories
which are informative of the agent’s behavior under a test-time state
distribution. To achieve this, we first define a prior over the types of
states we expect to see at test time. Unfortunately, most simple priors (e.g.
uniform) will also include many unreachable states which are not informative
to show the user. To avoid this issue, we use an exploration policy to seek
out states which match our prior distribution. By navigating to these new
states rather than directly initializing our environment in these states, we
reach a diverse yet still reachable set of states.

To summarize our contributions, we designed counterfactual trajectories that
explain behavior policies in out-of-distribution states and conducted studies
with human users that show increased understanding and generalization to out-
of-distribution states by measuring prediction ability. 111Our code is
available at https://github.com/juliusfrost/cfrl-rllib

Figure 2: Explanation pipeline

##  2 Related Work

###  2.1 Saliency

One popular class of agent interpretability approaches is saliency methods,
which aim to show which features from the input cause specific agent output.
This method was employed by (Greydanus et al., 2017), (Hilton et al., 2020),
and (Puri et al., 2020) to interpret which parts of the input that were deemed
important to the agent’s decision. (Anderson et al., 2019) use saliency as a
part of a more general agent explanation method. Even with all of the success,
saliency methods have been shown to have limitations. (Atrey et al., 2020)
show that saliency does not necessarily correspond to underlying
representation in RL. In addition, saliency methods are not designed to deal
with agents with memory, and in general may not be informative in a multi-
timestep environment. All of these challenges point to saliency as a helpful
explanatory tool but not as a stand-alone solution.

###  2.2 Interpretable Representations

Another approach is to use interpretable intermediate representations. In this
approach, models are constructed so that even if they are not interpretable
end-to-end, there will be some interpretable bottleneck or set of features
used to make the final decision. (Madumal et al., 2019) use a causal model to
generate explanations. (Chen et al., 2015) apply this method to autonomous
driving, where a model is trained to take an image and output understandable
intermediate representations in addition to a final action. Other work such as
(Kim et al., 2018) and (Jiang et al., 2019) use language as an informative
intermediate representation or output. While successful in some applications,
these kinds of interpretability methods often rely on domain knowledge to
choose appropriate intermediates, and they are not applicable if a policy’s
decision cannot be explained by a small set of interpretable components.

###  2.3 Critical States

In the critical states framework, explanations are generated by showing
especially informative “critical” states or trajectories from an agent. For
example the work of (Huang et al., 2018) considers states critical if they
have a large Q-value-difference between actions or a very low entropy action
distribution under a maximum entropy learning regimen. (Amir & Amir, 2018)
employ a similar approach, but they use videos instead of states and filter
for diversity of trajectories. These methods can be seen as extensions over a
baseline method of visualizing random states or failure states of an agent in
an environment. The explanation format of our method - videos of the agent’s
behavior in informative situations - is closest to this line of work.

###  2.4 Counterfactual Style Explanations

In (Rupprecht et al., 2019), researchers build a generative model of states
and try to show users the agent’s behavior in states that optimize some notion
of "interesting." Again using a generative model, (Olson et al., 2019) use a
GAN to produce modifications to states such that in the new state an agent
takes a different action than it would have in the original state. While these
works synthesize states directly, our method uses valid actions to access
different parts of the state space. This means that any state we visit is
actually reachable, at least in the training environment. The closest method
to ours appears in the work of (Witty et al., 2018), in which the authors try
to characterize generalization by exploring different starting configurations,
either by directly changing the start state or by using states visited by
another agent. In contrast to our work, their insights are applied to
characterizing generalization, and not as explanations for humans to
understand agent behavior.

(Anderson et al., 2019) provide a detailed empirical user study on the
effectiveness of saliency map and reward explanations. They conclude there is
no one explanation that fits all instances, but using several methods yields
the best mental model, and thus encourage using multiple methods for
explainability. We hope to introduce a novel explainability method which may
add to this suite of tools.

##  3 Method

###  3.1 Preliminaries

In this work, we consider the challenge of generating informative trajectories
which help a human user understand an agent’s policy (called the behavioral
policy) $\pi_{\theta}(a_{t}|s_{t})$, which is parametrized by weights $\theta$
which can be trained through any RL or imitation learning method. The agent
takes actions $a_{t}$ in states $s_{t}$, resulting in a trajectory $\tau$.

Figure 3: Distribution shifts, illustrated in the MiniGrid environment. Our
user study was performed in an environment similar to this, except with the
bottom-right room connected to the others.

Our goal is to help the user predict the agent’s behavior in a new test
environment in the presence of a state distribution shift which results in
different trajectory probabilities at train and test time:
$p^{\text{train}}(\tau|\theta)\neq p^{\text{test}}(\tau|\theta)$. We consider
two types of environment changes which can produce this distribution shift:
differences in the initial state distribution
($p_{0}^{\text{train}}(s_{0})\neq p_{0}^{\text{test}}(s_{0})$ for some initial
state $s_{0}$) and slight differences in the dynamics
($p^{\text{train}}(s_{t}|s_{t-1},a_{t-1})\neq
p^{\text{test}}(s_{t}|s_{t-1},a_{t-1})$ for some state $s_{t}$ and action
$a_{t}$). We assume that there are no changes to the agent’s action space or
policy at test time - i.e. $\pi_{\theta}(a_{t}|s_{t})$ remains consistent
between train and test time. We also assume that dynamics distribution shifts
are minor. See Section 3.4 for more detailed discussion on this.

###  3.2 Method Overview

Figure 2 summarizes our method. We build upon the explanatory pipeline used in
past critical states work such as (Huang et al., 2018). An agent is first
pretrained using any RL or imitation learning algorithm. We then collect many
trajectory rollouts of the agent, select a set to show the user, and finally
test whether these explanatory trajectories have helped the human better
predict the agent’s behavior in the test environment. This is challenging
because of the distribution shift between the train and test environments.

Past work has attempted to select an informative set of trajectories to show
the user by innovating on Step 3 in Figure 2 \- the trajectory selection
process. The critical states line of work selects trajectory segments which
show the agent’s behavior at states where the policy believes one particular
action is much better than the others (Huang et al., 2018). Unfortunately,
methods like these which sub-select trajectories from a set of agent rollouts
in the train environment are often not very informative about the agent’s
behavior under a distribution shift. Consider the train/test initial state
distribution shift in Figure 3. A well-trained agent which was always
initialized in the top-left room during training may take a nearly-
deterministic trajectory navigating toward the goal every time, so none of the
trajectories visited by the agent in the training environment will ever visit
the states in the top-right room which the agent will experience at test time.

One way to generate trajectories in the training distribution which are more
informative of the test distribution would be to manually change the starting
state distribution. For instance, in the grid world shown in Figure 3, the
simulator could be modified to directly initialize the agent in a uniform
distribution over starting cells, ensuring that some trajectories will have
been seen in whatever state the agent finds itself at test time.

One challenge with this approach is that it is often difficult to restrict the
set of randomized initialization to reachable states the agent could possibly
visit at test time. Instead, we propose creating a new distribution of start
states using an exploratory policy $\pi_{\phi}$ which can navigate from the
states seen on the training distribution to a new, more diverse distribution
of start states $p_{0}^{\text{expl}}(s_{0}|\phi)$ which are guaranteed to be
reachable. Implementing this only involves adding one additional Step 1.5 to
the pipeline shown in Figure 2.

Our trajectory generation method is summarized in the following algorithm:

  1. 1.

Run the behavioral policy $\pi_{\theta}$ in the train environment. Pause at a
randomly selected state along the trajectory.

  2. 2.

Run exploratory policy $\pi_{\phi}$ starting at this state for a fixed number
of timesteps. The exploratory policy will guide the agent to a set of states
not typically observed in the training data. We call these counterfactual
states because we can frame it as a counterfactual question: “what would the
behavioral policy do if it instead navigated to this state?”

  3. 3.

Run the behavioral policy starting from the counterfactual state for the
remainder of the trajectory.

  4. 4.

Show the user video rollouts of the agent’s behavior starting in these
counterfactual states.

Figure 4: Possible types of train/test distribution shifts. Our user
experiments use Shift A, although we would expect our method to be applicable
to Shifts C and E as well.

Note that the trajectories selection method used in item 4 of our algorithm
can be very simple - in our experiments we randomly select trajectories.
However, our method is complementary to other work on selecting informative
trajectories, so future work could select trajectories more systematically,
for instance filtering out near-repeated trajectories as was done in (Amir &
Amir, 2018) or selecting based on critical states as was done in (Huang et
al., 2018).

###  3.3 Exploration Objective

Our aim is to choose an exploration objective which will guide the agent to
states which produce maximally informative trajectories. Intuitively,
explanatory trajectories will be more informative the closer the distribution
of trajectories visited by our exploratory policy is to the test-time
distribution of trajectories visited. More formally, let
$p_{0}^{\text{expl}}(s_{0}|\phi)$ be the distribution of states in which the
agent finds itself after running the exploration policy, and let
$p^{\text{expl}}(\tau|\phi,\theta)$ be the distribution of trajectories
generated by rolling out the behavioral policy $\pi_{\theta}\ss$ initialized
in $p_{0}^{\text{expl}}(s_{0}|\phi)$. Our goal is to achieve
$p^{\text{expl}}(\tau|\phi,\theta)\approx p^{\text{test}}(\tau|\theta)$.

It is easiest to choose a principled exploration objective by making the
assumption that train and test-time dynamics are approximately identical. In
this situation, the train-test distribution shift consists entirely of changes
to the initial state: $p_{0}^{\text{train}}(s_{0})\neq
p_{0}^{\text{test}}(s_{0})$. However, if we are able to choose an exploration
policy which seeks out states matching the test time distribution - i.e. if
$p_{0}^{\text{expl}}(s_{0}|\phi)\approx p_{0}^{\text{test}}(s_{0})$ \- then
because we roll out the same behavioral policy in environments with identical
dynamics starting from nearly identical start distributions, we will achieve
$p^{\text{expl}}(\tau|\phi,\theta)\approx p^{\text{test}}(\tau|\theta)$.

We do not precisely know $p_{0}^{\text{test}}(s_{0})$ before the agent
experiences the test environment, so instead we choose a prior for our test
time start state distribution. In the absence of domain knowledge about what
states are likely under the test-time distribution shift, we assume a uniform
prior over reachable test-time states and try to match this distribution as
closely as possible by using an exploratory policy which seeks out a uniform
distribution of states. (Note that with this exploration objective, the only
reason we begin to run $\pi_{\phi}$ along a trajectory generated by
$\pi_{\theta}$ rather than than running $\pi_{\phi}$ from the start is to make
achieving the exploration objective easier for a learned exploration policy.)
However, if the user does have domain knowledge about which states are likely
at test time, the exploratory policy could be modified to preferentially seek
out those states. For instance, if the user believes that test-time states
will be near the training distribution, the exploration objective may instead
be to seek out a uniform distribution of states within a certain radius of
those seen at train time.

There has been a variety of past exploration work in which agents are trained
to seek out a uniform distribution of states or to learn distinct skills which
take the agent to different parts of the state space such as (Sharma et al.,
2019) and (Eysenbach et al., 2018). While any of these off-the-shelf
exploration algorithms could be used, in our work, we experiment in the
Minigrid environment, where the state space and dynamics are simple enough it
is possible to hard-code an oracle exploration policy which achieves a uniform
coverage over reachable states. By using this oracle policy, we are able to
directly test the question of whether highly diverse states serve as better
explanations rather than implicitly also testing the ability of the
exploration policy.

###  3.4 Distribution shifts where Our method helps

There are multiple types of changes between train and test environments which
can induce a difference between $p^{\text{train}}(\tau|\theta)$ and
$p^{\text{test}}(\tau|\theta)$, as illustrated in Figure 4. In Section 3.3, we
gave an intuitive justification for why we would expect our method to help in
cases where train and test-time dynamics are identical, and where the
exploration policy is able to cover the test-time start distribution - Shift
A, in Figure 4. Now, we consider our method’s usefulness in the presence of
other distribution shifts.

Shift B starts the agent in a new state which is unreachable in the train
environment, so it is unlikely that our explanations can help there. This
restriction to only reachable distribution shifts is very limiting, but we can
partially relax it in cases where the distribution shift does not change the
agent’s internal representation of the environment. Shift C illustrates an
example of this: the test-time environment has the top-right corner wall
segment removed, but so long as the agent’s internal representation of the
environment does not change much based on this change, the agent will act the
same way it would act if the wall was present. As a result, explanations
generated in the train environment will still be informative about test-env
behavior.

Finally, we consider dynamics shifts. In the case of global dynamics shifts,
as shown in Shift D, our method is likely unhelpful, since train-time and
test-time trajectories look very different even when starting from the same
state. However, if dynamics shifts only occur in a few states in the test
environment as shown in Shift E, then most of the time a train-time trajectory
and a test-time trajectory beginning at the same state will look identical.
Therefore, we can think of the dynamics shift purely as a way to introduce a
state distribution shift. While our explanations will not help the user
predict these dynamics shifts, the user should be able to predict how the
agent will act after the dynamics shift leads the agent to an unseen state.

##  4 Experiments

In this section, we test our primary hypothesis that diverse states help the
user understand the behavioral policy’s performance in the presence of a
distribution shift. We empirically evaluate this through a Mechanical Turk
user study for which we obtained IRB approval.

###  4.1 Study Design

Our study consisted of a questionnaire hosted online. Our participants were
Mechanical Turkers from the United States with a high Mechanical Turk approval
rating on past jobs. Each questionnaire consisted of an explanation phase in
which the user sees video trajectories of the policy in the train environment
and an evaluation phase which measures the user’s understanding of the
policy’s performance in the test environment. Train and test environments
differ in which room the agent is initialized, as shown in Figure 5. We tested
participants on two tasks, described in Sections 4.1.2 and 4.1.3.

####  4.1.1 Explanation Types

For each of the tasks, our user study compares three different types of
explanations. In the _Random States_ setting, participants are shown 10 sample
trajectories selected at random from a dataset of rollouts in an environment
with the agent initialized in a particular start region. In the _Critical
States_ setting, participants are shown a trajectory containing each of the 10
lowest entropy states in the dataset, all with the agent initialized in the
same starting region each time. If two or more of the low-entropy states are
in the same trajectory, the next lowest entropy states are selected. Finally,
in the _Counterfactual States_ setting, participants are shown 10 trajectories
of the behavioral policy beginning from a state chosen by the oracle
exploration policy, which seeks out a uniform distribution over reachable
states. The _Random States_ setting acts as a baseline for the other two
methods, because it is common practice to look at random rollouts to
understand policy behavior. The _Critical States_ setting was chosen as
another baseline, because it is the prior work with an explanation format
which most closely matches our explanations, allowing for easier comparison.
Finally, the _Counterfactual States_ condition is our proposed explanation
method. Since the agent is initialized in a different distribution of start
states at test time, we expect the _Random States_ and _Critical States_
methods to be uninformative at helping the user predict test-time performance.
In contrast, _Counterfactual States_ show more diversity and should enable to
the user to predict test-time performance.

####  4.1.2 Behavior Understanding Task

In this task, users first see videos of rollouts of $\pi_{\theta}$ in the
training environment selected using one of the three explanation methods
described in section 4.1.1. In this task, as users watch these explanations,
they are asked to build a mental model of the specific behaviors of
$\pi_{\theta}$. They are then presented with a context state in the test
environment, which is identical to the train environment except that the agent
is initialized in a different part of the state space where $\pi_{\theta}$
performs poorly. Below the context video, the user is presented with videos of
three potential continuation trajectories from the context state. Only one of
the continuation trajectories is generated by $\pi_{\theta}$, and participants
must select the continuation corresponding to $\pi_{\theta}$. The incorrect
choices for continuation trajectories were generated with policies manually
designed to have behavior which is distinctly different from $\pi_{\theta}$ \-
one of which is always successful in the test distribution, and one which
fails in a way which is visually distinct from $\pi_{\theta}$.

####  4.1.3 Performance Evaluation Task

We next evaluate whether the user can predict the performance of
$\pi_{\theta}$ in the test environment. To do this, users are shown a context
state from a test environment where the agent starts in a different state. The
user must guess whether $\pi_{\theta}$ will succeed or fail from that state,
where success means the agent ends in the desired goal location. This task is
designed to measure the participants’ understanding of the probability an
agent succeeds in its task given its state.

####  4.1.4 Experiment Settings

We test our method in a custom Minigrid environment (Chevalier-Boisvert et
al., 2018) where there are four rooms with doors connecting between them
(shown in Figure 4A). We train our behavioral policy using the A2C algorithm
(Mnih et al., 2016) to good performance on the train environment. To select
counterfactual states, we use an oracle policy which uniformly samples a
feasible state and navigates towards it. In the study, each participant sees
10 explanations and then answers 10 questions of our chosen evaluation task.

For the Behavior Understanding task, the policy was trained in the top-left
room, and explanations were generated in this room. In the test environment,
the agent is positioned in the bottom-right room. As a result, the agent
typically succeeds when initialized in the top-left (where we collect
explanations) but typically fails when initialized in the top-right room
(where we initialize in the evaluation phase).

For the Policy Evaluation task, we test users on two different behavioral
policies. The first is only able to succeed starting in the top-left room, and
the other is only able to succeed from the bottom-right room. In both cases,
we collect explanations in the top-left room. Users are evaluated on the
agent’s behavior in all four rooms.

Figure 5: Experimental setup for the counterfactual states user study. Left: users see 10 train-time rollouts such as these in the explanation phase; Right: users are asked to either select the test-time behavior which corresponds with the policy they saw previously (Task 1 - Behavior Understanding) or predict the agent’s success in a new state (Task 2 - Policy Evaluation). Table 1: Minigrid Policy Understanding and Evaluation | Task 1 |  | Task 2 |
---|---|---|---|---
Condition | Accuracy | n | Accuracy | n
Random | 0.04 | 10 | 0.6315 | 19
Critical | 0.1458 | 10 | 0.6411 | 17
Counterfactual | 0.3638 | 10 | 0.6904 | 21
Table 2: Task 1 is the Behavior Understanding task. Task 2 is Performance
Evaluation. Bold means the result is statistically greater than the rest for
$p<0.05$ with a one-sided T-test.

###  4.2 Results

Qualitatively, in both tasks, the counterfactual states method produces more
informative rollouts. In both tasks, neither the random states nor critical
states methods produce any explanatory trajectories where the agent visits the
room where the it is initialized at test time. This makes it hard for the user
to predict the agent’s performance in these states. However, in the
counterfactual states condition, the user sees multiple examples of the agent
in the distribution of states it will see at test time.

Quantitative results, however, are mixed. Table 2 shows that for the Behavior
Understanding task, users are able to select the correct behavior $36.38\%$ of
the time, a statistically significant improvement over the random and critical
states conditions. This performance is only marginally above random guessing
($33.3\%$), but at least the user did not seem to be actively misled by the
explanations, as occurred in the random and critical states cases. On the
Policy Evaluation task, users perform similarly in all conditions, with no
statistically significant difference between them.

##  5 Discussion

While our method outperformed baselines in one of two tasks, there is
significant room for improvement, as users did not consistently perform well
in any condition of any task, even though the environment and agent behavior
were both quite simple. Part of the difficulty in constructing good
explanations derives from the population we were testing. Mechanical Turkers
likely have little prior experience with reinforcement learning or time to
carefully analyze agent behavior. Anecdotally, we also observed that many
survey-takers come into the task with their own biases and preconceptions (for
instance, assuming that a policy will be deterministic rather than stochastic,
or assuming that if an agent does well in one situation it probably does well
everywhere). Future work could test the usefulness of our counterfactual
states method for agent developers who can spend time to familiarize
themselves with the task, the explanation interface, and the types of
behaviors which typically emerge in policies. In this setting, we could also
test our method’s usefulness in more complex environments. The experiments we
have run so far are quite limited: the environment is simple, an oracle
exploration policy is available, the distribution shift only changes the
agent’s starting position, and it is easy for a user to distinguish between
correct and incorrect behaviors. Future work could explore whether
counterfactual states hold value in more realistic settings.

In this paper, we hypothesized that showing higher-diversity reachable
trajectories will enable users to better understand agent behavior. We tested
this hypothesis through a user study in a gridworld in which the test-time
initial state distribution differed from the train-time start distribution.
Our method showed improvements over less diverse baseline methods in one of
two tasks, but our experiments show room for improvement and illustrate the
ongoing challenge of creating explanations which improve the understanding of
non-expert users.

##  6 Acknowledgements

Thanks to Anna Rohrbach for ongoing feedback on this work.

## References

  * Akkaya et al. (2019) Akkaya, I., Andrychowicz, M., Chociej, M., Litwin, M., McGrew, B., Petron, A., Paino, A., Plappert, M., Powell, G., Ribas, R., et al.  Solving rubik’s cube with a robot hand.  _arXiv preprint arXiv:1910.07113_ , 2019.
  * Amir & Amir (2018) Amir, D. and Amir, O.  Highlights: Summarizing agent behavior to people.  In _Proceedings of the 17th International Conference on Autonomous Agents and MultiAgent Systems_ , pp. 1168–1176, 2018.
  * Anderson et al. (2019) Anderson, A., Dodge, J., Sadarangani, A., Juozapaitis, Z., Newman, E., Irvine, J., Chattopadhyay, S., Fern, A., and Burnett, M.  Explaining reinforcement learning to mere mortals: An empirical study.  _IJCAI_ , 2019.
  * Atrey et al. (2020) Atrey, A., Clary, K., and Jensen, D.  Exploratory not explanatory: Counterfactual analysis of saliency maps for deep rl.  _ICLR_ , 2020.
  * Bellemare et al. (2013) Bellemare, M. G., Naddaf, Y., Veness, J., and Bowling, M.  The arcade learning environment: An evaluation platform for general agents.  _Journal of Artificial Intelligence Research_ , 47:253–279, jun 2013.
  * Chen et al. (2015) Chen, C., Seff, A., Kornhauser, A., and Xiao, J.  Deepdriving: Learning affordance for direct perception in autonomous driving.  In _Proceedings of the IEEE international conference on computer vision_ , pp. 2722–2730, 2015.
  * Chevalier-Boisvert et al. (2018) Chevalier-Boisvert, M., Willems, L., and Pal, S.  Minimalistic gridworld environment for openai gym.  https://github.com/maximecb/gym-minigrid, 2018.
  * Eysenbach et al. (2018) Eysenbach, B., Gupta, A., Ibarz, J., and Levine, S.  Diversity is all you need: Learning skills without a reward function.  _arXiv preprint arXiv:1802.06070_ , 2018.
  * Eysenbach et al. (2020) Eysenbach, B., Asawa, S., Chaudhari, S., Levine, S., and Salakhutdinov, R.  Off-dynamics reinforcement learning: Training for transfer with domain classifiers.  2020\.
  * Greydanus et al. (2017) Greydanus, S., Koul, A., Dodge, J., and Fern, A.  Visualizing and understanding atari agents.  _ICML_ , 2017.
  * Hilton et al. (2020) Hilton, J., Cammarata, N., Carter, S., Goh, G., and Olah, C.  Understanding rl vision.  _Distill_ , 2020.  doi: 10.23915/distill.00029.  https://distill.pub/2020/understanding-rl-vision.
  * Hoffman et al. (2018) Hoffman, R. R., Mueller, S. T., Klein, G., and Litman, J.  Metrics for explainable ai: Challenges and prospects.  _arXiv preprint arXiv:1812.04608_ , 2018.
  * Huang et al. (2018) Huang, S. H., Bhatia, K., Abbeel, P., and Dragan, A. D.  Establishing appropriate trust via critical states.  In _2018 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_ , pp. 3929–3936, Oct 2018.  doi: 10.1109/IROS.2018.8593649.
  * Jiang et al. (2019) Jiang, Y., Gu, S., Murphy, K., and Finn, C.  Language as an abstraction for hierarchical deep reinforcement learning.  _arXiv preprint arXiv:1906.07343_ , 2019.
  * Kim et al. (2018) Kim, J., Rohrbach, A., Darrell, T., Canny, J., and Akata, Z.  Textual explanations for self-driving vehicles.  In _Proceedings of the European conference on computer vision (ECCV)_ , pp. 563–578, 2018.
  * Madumal et al. (2019) Madumal, P., Miller, T., Sonenberg, L., and Vetere, F.  Explainable reinforcement learning through a causal lens.  _AAAI 2020_ , 2019.
  * Mnih et al. (2016) Mnih, V., Badia, A. P., Mirza, M., Graves, A., Lillicrap, T., Harley, T., Silver, D., and Kavukcuoglu, K.  Asynchronous methods for deep reinforcement learning.  In _International conference on machine learning_ , pp. 1928–1937. PMLR, 2016.
  * Olson et al. (2019) Olson, M. L., Neal, L., Li, F., and Wong, W.-K.  Counterfactual states for atari agents via generative deep learning.  _IJCAI_ , 2019.
  * Peng et al. (2018) Peng, X. B., Andrychowicz, M., Zaremba, W., and Abbeel, P.  Sim-to-real transfer of robotic control with dynamics randomization.  In _2018 IEEE international conference on robotics and automation (ICRA)_ , pp. 3803–3810. IEEE, 2018.
  * Puri et al. (2020) Puri, N., Verma, S., Gupta, P., Kayastha, D., Deshmukh, S., Krishnamurthy, B., and Singh, S.  Explain your move: Understanding agent actions using specific and relevant feature attribution.  In _ICLR_ , 2020.
  * Rupprecht et al. (2019) Rupprecht, C., Ibrahim, C., and Pal, C. J.  Finding and visualizing weaknesses of deep reinforcement learning agents.  _arXiv preprint arXiv:1904.01318_ , 2019.
  * Rusu et al. (2017) Rusu, A. A., Večerík, M., Rothörl, T., Heess, N., Pascanu, R., and Hadsell, R.  Sim-to-real robot learning from pixels with progressive nets.  In _Conference on Robot Learning_ , pp. 262–270. PMLR, 2017.
  * Sharma et al. (2019) Sharma, A., Gu, S., Levine, S., Kumar, V., and Hausman, K.  Dynamics-aware unsupervised discovery of skills.  _arXiv preprint arXiv:1907.01657_ , 2019.
  * Tassa et al. (2018) Tassa, Y., Doron, Y., Muldal, A., Erez, T., Li, Y., de Las Casas, D., Budden, D., Abdolmaleki, A., Merel, J., Lefrancq, A., Lillicrap, T., and Riedmiller, M.  DeepMind control suite.  Technical report, DeepMind, January 2018.  URL https://arxiv.org/abs/1801.00690.
  * Witty et al. (2018) Witty, S., Lee, J. K., Tosch, E., Atrey, A., Littman, M., and Jensen, D.  Measuring and characterizing generalization in deep reinforcement learning.  _arXiv preprint arXiv:1812.02868_ , 2018.


&+\int\frac{d^{2}\mathbf{q}}{(2\pi)^{2}}\frac{d^{2}\mathbf{p}}{(2\pi)^{2}}\frac{\mathbf{q}\times\mathbf{k}}{|\mathbf{k}-\mathbf{q}|^{2}q_{\perp}^{2}}\frac{\mathbf{q}\cdot(\mathbf{q}-2\mathbf{p})}{4p_{\perp}^{2}|\mathbf{q}-\mathbf{p}|^{2}}\int_{-\pi}^{\pi}\frac{d\phi}{2\pi}\frac{1}{w_{\perp}^{2}}\Big{(}2\mathbf{p}\cdot(\mathbf{q}-\mathbf{p})\left[[b_{\perp}(\mathbf{p}),b_{\perp}(\mathbf{q}-\mathbf{p})],b_{\perp}(\mathbf{k}-\mathbf{q})\right]\\\
&\qquad\qquad-(q_{\perp}^{2}-w_{\perp}^{2}-2\mathbf{p}\cdot(\mathbf{q}-\mathbf{p}))\left[[b_{\eta}(\mathbf{p}),b_{\eta}(\mathbf{q}-\mathbf{p})],b_{\perp}(\mathbf{k}-\mathbf{q})\right]\Big{)}\\\
&\qquad\qquad\times\left\\{(w_{\perp}^{2}-|\mathbf{k}-\mathbf{q}|^{2})\int_{-\pi}^{\pi}\frac{d\theta}{2\pi}\frac{1}{u^{2}_{\perp}}(1-J_{0}(u_{\perp}\tau))+(1-J_{0}(|\mathbf{k}-\mathbf{q}|\tau))\right\\}\\\
&+\int\frac{d^{2}\mathbf{q}}{(2\pi)^{2}}\frac{d^{2}\mathbf{p}}{(2\pi)}\frac{\mathbf{q}\times\mathbf{p}}{p_{\perp}^{2}|\mathbf{q}-\mathbf{p}|^{2}}\Big{[}[b_{\perp}(\mathbf{p}),b_{\eta}(\mathbf{q}-\mathbf{p})],b_{\eta}(\mathbf{k}-\mathbf{q})\Big{]}\int_{-\pi}^{\pi}\frac{d\phi}{2\pi}\frac{\mathbf{q}\cdot(\mathbf{q}-2\mathbf{p})+w_{\perp}^{2}}{q_{\perp}^{2}-w_{\perp}^{2}}\\\
&\quad\times\left\\{\frac{w_{\perp}^{2}-|\mathbf{k}-\mathbf{q}|^{2}}{w_{\perp}|\mathbf{k}-\mathbf{q}|}\int_{-\pi}^{\pi}\frac{d\theta}{2\pi}\frac{\cos\theta}{u^{2}_{\perp}}(1-J_{0}(u_{\perp}\tau))-\frac{q_{\perp}^{2}-|\mathbf{k}-\mathbf{q}|^{2}}{q_{\perp}|\mathbf{k}-\mathbf{q}|}\int_{-\pi}^{\pi}\frac{d\phi^{\prime}}{2\pi}\frac{\cos\phi^{\prime}}{w_{\perp}^{\prime
2}}(1-J_{0}(w^{\prime}_{\perp}\tau))\right\\}\,.\\\ \end{split}$ (134)

Let us summarize and comment on the general procedures for solving the
classical Yang-Mills equations perturbatively in the dilute-dense regime. At
each fixed order $g^{2m+1}$, the dynamical fields
$\beta^{(2m+1)}(\tau,\mathbf{k})$ and
$\beta^{(2m+1)}_{\perp}(\tau,\mathbf{k})$ satisfy the inhomogeneous Bessel
differential equations of orders one and zero, respectively. They are solved
using the well-established method of variation of parameters. The success of
this method relies on recombining the time-dependent factors in the source
terms $S_{\eta}^{(2m+1)}$ and $S_{\perp}^{(2m+1)}$ so that $S_{\eta}^{(2m+1)}$
only depends on $J_{1}(c\tau)$ and $S_{\perp}^{(2m+1)}$ only depends on
$J_{0}(c\tau)$. Owing to the Graf’ s formula, this is always possible. As for
$\Lambda^{(2m+1)}(\tau,\mathbf{k})$, it satisfies first order differential
equation with source term $S_{\Lambda}^{(2m+1)}$. All one needs to do is
express the time-dependent factors in $S_{\Lambda}^{(2m+1)}(\tau,\mathbf{k})$
in terms of Bessel functions $J_{1}(c\tau)$. Each time Graf’s formula is used,
an extra angular integral is introduced. Unfortunately, the number of terms in
the solutions increases dramatically as one goes to higher and higher
perturbative orders and the problem quickly becomes unmanageable analytically.

##  8 Discussions and Outlooks

In this paper, we have presented the first step towards completing the
calculations of the first saturation corrections to physical observables in
high energy proton-nucleus collisions. We solved the classical Yang-Mills
equations in the dilute-dense regime beyond the leading order. We explicitly
constructed the order-$g^{3}$ and order-$g^{5}$ solutions. The main results
are presented in Eqs. (70), (83), (95), (117), (126) and (134). The major
mathematical technique that makes the analytic solutions possible is Graf’s
formula, which expresses product of two Bessel functions in terms of an
angular integral of one Bessel function. As a consistence check of our main
results, when the target Wilson line $U(\mathbf{x})=1$, the gluon fields
vanish $\beta(\tau,\mathbf{k})=0$ and $\beta_{i}(\tau,\mathbf{k})=0$, i.e.
there is no gluon production in the absence of scattering as expected.

There are a few apparent features of the solutions. First of all, the time
dependent factors in the longitudinal gluon field

$\beta(\tau,\mathbf{k})=\beta^{(1)}(\tau,\mathbf{k})+\beta^{(3)}(\tau,\mathbf{k})+\beta^{(5)}(\tau,\mathbf{k})+\ldots$
(135)

are uniquely determined by one single type of Bessel functions $J_{1}(a\tau)$
although the values of the argument $a$ might be different at different
orders. On the other hand, the time dependent factors in the transverse field

$\beta_{i}(\tau,\mathbf{k})=\beta^{(1)}_{i}(\tau,\mathbf{k})+\beta^{(3)}_{i}(\tau,\mathbf{k})+\beta^{(5)}_{i}(\tau,\mathbf{k})+\ldots$
(136)

are completely determined by $J_{0}(a\tau)$ again with possible different
arguments $a$.

Second, the order-$g$ solutions
$\beta^{(1)}(\tau,\mathbf{k}),\beta^{(1)}_{i}(\tau,\mathbf{k})$ do not involve
mutual interactions of the “valence” color charges in the proton. Each color
charge in the proton independently scatter on the target. On the other hand,
the order-$g^{3}$ solutions
$\beta^{(3)}(\tau,\mathbf{k}),\beta^{(3)}_{i}(\tau,\mathbf{k})$ represent
interactions of two color charges in the proton. Their interactions can happen
before or after the collisions with the target. The order-$g^{5}$ solutions
$\beta^{(5)}(\tau,\mathbf{k}),\beta^{(5)}_{i}(\tau,\mathbf{k})$ represent
interactions of three color charges in the proton. Their interactions likewise
can happen before or after the collisions with the target. It can also be that
two color charges interact with each other before scattering on the target and
then interact with the third color charge only after the collisions with the
target.

Another important feature is related to the gauge dependence of the concepts
of initial state effects and final state effects as defined on the basis of
the field. In the main context of the paper, the solutions are given in the
initial time Coulomb sub gauge. In appendix D, we present the order-$g$ and
order-$g^{3}$ solutions in the non-Coulomb sub gauge. These two sets of
solutions are related by a gauge transformation. Performing direct comparison
of these two sets of solutions, would convince the reader that some final
state effects in the non-Coulomb subgauge become the initial state effects in
the Coulomb subgauge.

On the other hand, physical observables are independent of gauge
transformations. With the solutions of the classical Yang-Mills equations at
hand, one can calculate several interesting physical quantities that can be
constructed from the classical gluon fields. For example, the the energy-
momentum tensor of the gluon fields produced in high energy pA collisions by

$T^{\mu\nu}(\tau,\mathbf{x})=F^{\mu\lambda}(\tau,\mathbf{x})F_{\lambda}^{\,\,\,\,\nu}(\tau,\mathbf{x})+\frac{1}{4}g^{\mu\nu}F^{\kappa\lambda}(\tau,\mathbf{x})F_{\kappa\lambda}(\tau,\mathbf{x}).$
(137)

Tracing over the color matrix is understood in the above definition like
$AB=A^{a}B^{a}=\mathrm{Tr}(AB)$ with $a=1,\ldots N_{c}^{2}-1$.

In the second paper of this series, we will calculate the single inclusive
soft gluon production and double inclusive soft gluon production. These
observables are constructed using the LSZ formula

$\begin{split}&\hat{a}_{\mathbf{p}}^{a\dagger}(\tau)=-i\tau\sqrt{\frac{\pi}{4}}\left(H_{1}^{(2)}(p_{\perp}\tau)\overleftrightarrow{\partial_{\tau}}\tilde{\beta}^{a}(\tau,\mathbf{p})\right),\\\
&\hat{c}_{\mathbf{p}}^{a\dagger}(\tau)=-i\tau\sqrt{\frac{\pi}{4}}\left(H_{0}^{(2)}(p_{\perp}\tau)\overleftrightarrow{\partial_{\tau}}\beta_{\perp}^{a}(\tau,\mathbf{p})\right),\end{split}$
(138)

and taking the limit $\tau\rightarrow\infty$. Here $H_{0}^{(2)}(x)$ and
$H_{1}^{(2)}(x)$ are Hankel functions of second kind, order zero and order
one, respectively. The left right derivative is defined as
$f_{1}(x)\overleftrightarrow{\partial_{x}}f_{2}(x)=f_{1}(x)\partial_{x}f_{2}(x)-\partial_{x}f_{1}(x)f_{2}(x)$.
Unlike the energy-momentum tensor which is gague invariant by construction,
the single inclusive gluon production by
$\hat{a}_{\mathbf{p}}^{a,\dagger}\hat{a}_{\mathbf{p}}^{a}+\hat{c}_{\mathbf{p}}^{a,\dagger}\hat{c}_{\mathbf{p}}^{a}$
are explicitly dependent on $\beta(\tau,\mathbf{k})$ and
$\beta_{\perp}(\tau,\mathbf{k})$ which are gauge variant objects. Thus gauge
invariance of the single inclusive production will serve as a non-trivial of
the derived solutions.

It should be mentioned that to solve the classical Yang-Mills equations in the
dense-dense regime, another semi-analytical approach was developed in Refs.
Fries:2006pv ; Chen:2015wia . In this method, the solutions are expanded as
power series expansions in proper time $\tau$. Although at each order in
$\tau$, all the saturation effects are included, the solutions are only
meaningful when the values of $\tau$ are small. This small-$\tau$ expansion
method cannot be used to rigorously calculate the single inclusive gluon
production by the LSZ formula which requires taking the
$\tau\rightarrow\infty$ limit. Our expansions in terms of coupling constant
$g$ are valid for all the proper time but can only take into account the
saturation effects order by order. Amusingly one can combine these two
analytical approaches into a double expansion in $\tau$ and $g$ which gives a
non-trivial insights into the dense-dense regime. We defer further discussion
for a future publication.

###### Acknowledgements.

We thank A. Dumitru, A. Kovner, Yu. Kovchegov, M. Lublinsky, L. McLerran, and
R. Venugopalan for insightful discussions and collaboration on related
projects. We thank H. Duan for his contribution in the exploratory stage of
this project. We acknowledge support by the DOE Office of Nuclear Physics
through Grant No. DE-SC0020081.

##  Appendix A The non-dynamical field $\Lambda(\tau,\mathbf{x})$

In this appendix, we give a general proof that $\Lambda(\tau,\mathbf{x})$ is
not a dynamical field. To be precise, we want to show that the second order
differential equation for $\Lambda(\tau,\mathbf{x})$ can be obtained from the
first order constraint equation. The proof is done perturbatively using
induction. The perturbative expansions for the solutions are

$\begin{split}&\tilde{\beta}(\tau,\mathbf{x})=\sum_{n=0}^{\infty}g^{2n+1}\tilde{\beta}^{(2n+1)}(\tau,\mathbf{x}),\\\
&\tilde{\beta}_{i}(\tau,\mathbf{x})=\sum_{n=0}^{\infty}g^{2n+1}\beta_{i}^{(2n+1)}(\tau,\mathbf{x}),\\\
&\Lambda(\tau,\mathbf{x})=\sum_{n=0}^{\infty}g^{2n+1}\Lambda^{(2n+1)}(\tau,\mathbf{x}).\\\
\end{split}$ (139)

Substituting these expansions into the classical Yang-Mills equations in Eq.
(42), at order-$g^{2N+1}$, the first order constraint equation is

$\partial_{\tau}\partial^{2}\Lambda^{(2N+1)}=i\sum_{M=0}^{N-1}[\beta_{i}^{(2N-2M-1)},\partial_{\tau}\beta_{i}^{(2M+1)}]+i\sum_{M=0}^{N-1}[\tilde{\beta}^{(2N-2M-1)},\partial_{\tau}\tilde{\beta}^{(2M+1)}]\,.$
(140)

The right hand side of this equation involves all the lower order solutions
$\beta_{i}^{(K)},\tilde{\beta}^{(K)}$ for $K=1,3,\dots,2N-1$.

The second order differential equation for the order-$g^{2N+1}$ field is

$\begin{split}&\tau^{2}\partial_{\tau}^{2}\partial^{2}\Lambda^{(2N+1)}+\tau\partial_{\tau}\partial^{2}\Lambda^{(2N+1)}\\\
=&i\tau^{2}\sum_{M=0}^{N-1}[\tilde{\beta}^{(2N-2M-1)},\partial^{2}\tilde{\beta}^{(2M+1)}]+i\tau^{2}\sum_{M=0}^{N-1}[\beta_{i}^{(2N-2M-1)},(\partial^{2}\delta_{ij}-\partial_{i}\partial_{j})\beta_{j}^{(2M+1)}]\\\
&+\tau^{2}\sum_{M=0}^{N-1}\sum_{L=0}^{M-1}\partial_{i}[\tilde{\beta}^{(2N-2M-1)},[\beta_{i}^{(2M-2L-1)},\tilde{\beta}^{(2L+1)}]\\\
&+\tau^{2}\sum_{M=0}^{N-1}\sum_{L=0}^{M-1}\partial_{i}[\beta_{j}^{(2N-2M-1)},[\beta_{i}^{(2M-2L-1)},\beta_{j}^{(2L+1)}]].\end{split}$
(141)

Our goal is to prove that from the order-$g^{2N+1}$ constraint equation Eq.
(140), using all the lower order equations of motion, the second order
differential equation Eq. (141) can be derived. Taking time derivative of Eq.
(140) w.r.t. $\tau$, one obtains

$\tau^{2}\partial^{2}_{\tau}\partial^{2}\Lambda^{(2N+1)}=i\sum_{M=0}^{N-1}[\beta_{i}^{(2N-2M-1)},\tau^{2}\partial^{2}_{\tau}\beta_{i}^{(2M+1)}]+i\sum_{M=0}^{N-1}[\tilde{\beta}^{(2N-2M-1)},\tau^{2}\partial^{2}_{\tau}\tilde{\beta}^{(2M+1)}].$
(142)

Combining Eqs. (140) and (142), we obtain

$\begin{split}&\tau^{2}\partial_{\tau}^{2}\partial^{2}\Lambda^{(2N+1)}+\tau\partial_{\tau}\partial^{2}\Lambda^{(2N+1)}=i\sum_{M=0}^{N-1}[\beta_{i}^{(2N-2M-1)},\tau^{2}\partial^{2}_{\tau}\beta_{i}^{(2M+1)}+\tau\partial_{\tau}\beta_{i}^{(2M+1)}]\\\
&\qquad+i\sum_{M=0}^{N-1}[\tilde{\beta}^{(2N-2M-1)},\tau^{2}\partial^{2}_{\tau}\tilde{\beta}^{(2M+1)}+\tau\partial_{\tau}\tilde{\beta}^{(2M+1)}]\\\
\end{split}$ (143)

It contains lower order equations of motion. They are (for $M=0,1,\ldots,N-1$)

$\begin{split}&\tau^{2}\partial^{2}_{\tau}\tilde{\beta}^{(2M+1)}+\tau\partial_{\tau}\tilde{\beta}^{(2M+1)}=\tilde{\beta}^{(2M+1)}+\tau^{2}\partial_{i}^{2}\tilde{\beta}^{(2M+1)}\\\
&\qquad-i\tau^{2}\sum_{M^{\prime}=0}^{M-1}\partial_{i}[\beta_{i}^{(2M-2M^{\prime}-1)},\tilde{\beta}^{(2M^{\prime}+1)}]-i\tau^{2}\sum_{M^{\prime}=0}^{M-1}[\beta_{i}^{(2M-2M^{\prime}-1)},\partial_{i}\tilde{\beta}^{(2M^{\prime}+1)}]\\\
&\qquad-\tau^{2}\sum_{M^{\prime}=0}^{M-1}\sum_{N^{\prime}=0}^{M^{\prime}-1}[\beta_{i}^{(2M-2M^{\prime}-1)},[\beta_{i}^{2M^{\prime}-2N^{\prime}-1},\tilde{\beta}^{2N^{\prime}+1}]]\\\
\end{split}$ (144)

and

$\begin{split}&\tau^{2}\partial^{2}_{\tau}\beta_{i}^{(2M+1)}+\tau\partial_{\tau}\beta_{i}^{(2M+1)}\\\
=&\tau^{2}(\partial^{2}\delta_{ij}-\partial_{i}\partial_{j})\beta_{j}^{(2M+1)}+i\tau^{2}\sum_{M^{\prime}=0}^{M-1}[\tilde{\beta}^{(2M-2M^{\prime}-1)},\partial_{i}\tilde{\beta}^{(2M^{\prime}+1)}]\\\
&-i\tau^{2}\sum_{M^{\prime}=0}^{M-1}\partial_{j}[\beta_{j}^{(2M-2M^{\prime}-1)},\beta_{i}^{(2M^{\prime}+1)}]-i\tau^{2}\sum_{M^{\prime}=0}^{M-1}[\beta_{j}^{(2M-2M^{\prime}-1)},\partial_{j}\beta_{i}^{(2M^{\prime}+1)}-\partial_{i}\beta_{j}^{(2M^{\prime}+1)}]\\\
&+\tau^{2}\sum_{M^{\prime}=0}^{M-1}\sum_{N^{\prime}=0}^{M^{\prime}-1}[\tilde{\beta}^{(2M-2M^{\prime}-1)},[\beta_{i}^{(2M^{\prime}-2N^{\prime}-1)},\tilde{\beta}^{(2N^{\prime}+1)}]]\\\
&+\tau^{2}\sum_{M^{\prime}=0}^{M-1}\sum_{N^{\prime}=0}^{M^{\prime}-1}[\beta_{j}^{(2M-2M^{\prime}-1)},[\beta_{i}^{(2M^{\prime}-2N^{\prime}-1)},\beta_{j}^{(2N^{\prime}+1)}]].\\\
\end{split}$ (145)

Substituting these two equations in Eq. (143), we get

$\begin{split}&\tau^{2}\partial_{\tau}^{2}\partial^{2}\Lambda^{(2N+1)}+\tau\partial_{\tau}\partial^{2}\Lambda^{(2N+1)}\\\
=&i\sum_{M=0}^{N-1}[\beta_{i}^{(2N-2M-1)},\tau^{2}\partial^{2}_{\tau}\beta_{i}^{(2M+1)}+\tau\partial_{\tau}\beta_{i}^{(2M+1)}]+i\sum_{M=0}^{N-1}[\tilde{\beta}^{(2N-2M-1)},\tau^{2}\partial^{2}_{\tau}\tilde{\beta}^{(2M+1)}+\tau\partial_{\tau}\tilde{\beta}^{(2M+1)}]\\\
=&i\sum_{M=0}^{N-1}[\beta_{i}^{(2N-2M-1)},\tau^{2}(\partial^{2}\delta_{ij}-\partial_{j}\partial_{j})\beta_{j}^{(2M+1)}]+i\sum_{M=0}^{N-1}[\tilde{\beta}^{(2N-2M-1)},\tau^{2}\partial_{i}^{2}\tilde{\beta}^{(2M+1)}]\\\
&-\tau^{2}\sum_{M=0}^{N-1}\sum_{M^{\prime}=0}^{M-1}[\beta_{i}^{(2N-2M-1)},[\tilde{\beta}^{(2M-2M^{\prime}-1)},\partial_{i}\tilde{\beta}^{(2M^{\prime}+1)}]]\\\
&+\tau^{2}\sum_{M=0}^{N-1}\sum_{M^{\prime}=0}^{M-1}[\beta_{i}^{(2N-2M-1)},\partial_{j}[\beta_{j}^{(2M-2M^{\prime}-1)},\beta_{i}^{(2M^{\prime}+1)}]]\\\
&+\tau^{2}\sum_{M=0}^{N-1}\sum_{M^{\prime}=0}^{M-1}[\beta_{i}^{(2N-2M-1)},[\beta_{j}^{(2M-2M^{\prime}-1)},\partial_{j}\beta_{i}^{(2M^{\prime}+1)}-\partial_{i}\beta_{j}^{(2M^{\prime}+1)}]]\\\
&+\tau^{2}\sum_{M=0}^{N-1}\sum_{M^{\prime}=0}^{M-1}[\tilde{\beta}^{(2N-2M-1)},\partial_{i}[\beta_{i}^{(2M-2M^{\prime}-1)},\tilde{\beta}^{(2M^{\prime}+1)}]]\\\
&+\tau^{2}\sum_{M=0}^{N-1}\sum_{M^{\prime}=0}^{M-1}[\tilde{\beta}^{(2N-2M-1)},[\beta_{i}^{(2M-2M^{\prime}-1)},\partial_{i}\tilde{\beta}^{(2M^{\prime}+1)}]]\\\
\end{split}$ (146)

First of all, we have used the fact that the following terms vanish

$\sum_{M=0}^{N-1}[\tilde{\beta}^{(2N-2M-1)},\tilde{\beta}^{(2M+1)}]=([\tilde{\beta},\tilde{\beta}])^{(2N)}=0$
(147)
$\displaystyle\sum_{M=0}^{N-1}\sum_{M^{\prime}=0}^{M-1}\sum_{N^{\prime}=0}^{M^{\prime}-1}[\beta_{i}^{(2N-2M-1)},[\beta_{j}^{(2M-2M^{\prime}-1)},[\beta_{i}^{(2M^{\prime}-2N^{\prime}-1)},\beta_{j}^{(2N^{\prime}+1)}]]]$
$\displaystyle=([\beta_{i},[\beta_{j},[\beta_{i},\beta_{j}]]])^{(2N-2)}=0$
(148)
$\begin{split}&\sum_{M=0}^{N-1}\sum_{M^{\prime}=0}^{M-1}\sum_{N^{\prime}=0}^{M^{\prime}-1}[\beta_{i}^{(2N-2M-1)},[\tilde{\beta}^{(2M-2M^{\prime}-1)},[\beta_{i}^{(2M^{\prime}-2N^{\prime}-1)},\tilde{\beta}^{(2N^{\prime}+1)}]]]\\\
&-\sum_{M=0}^{N-1}\sum_{M^{\prime}=0}^{M-1}\sum_{N^{\prime}=0}^{M^{\prime}-1}[\tilde{\beta}^{(2N-2M-1)},[\beta_{i}^{(2M-2M^{\prime}-1)},[\beta_{i}^{2M^{\prime}-2N^{\prime}-1},\tilde{\beta}^{2N^{\prime}+1}]]]\\\
=&([\beta_{i},[\tilde{\beta},[\beta_{i},\tilde{\beta}]]]-[\tilde{\beta},[\beta_{i},[\beta_{i},\tilde{\beta}]]])^{(2N-2)}=0\end{split}$
(149)

Trivial algebra shows that Eq. (146) reproduces Eq. (141). Q.E.D.

##  Appendix B The method of variation of parameters

For the sake of completeness, we briefly review the method of solving
inhomogeneous differential equations used in the main body of the paper. The
method is often called variation of parameters.

For a general second order inhomogeneous differential equation

$a(x)y^{\prime\prime}+b(x)y^{\prime}+c(x)y=f(x),$ (150)

there are four steps to obtain its general solutions:

  * •

find two independent solutions $y_{1}(x)$ and $y_{2}(x)$ of the homogeneous
equation

$a(x)y^{\prime\prime}+b(x)y^{\prime}+c(x)y=0\,.$ (151)

The general solution for the homogeneous equation is

$y_{h}(x)=C_{1}y_{1}(x)+C_{2}y_{2}(x)\,,$ (152)

where $C_{1}$ snd $C_{2}$ are coefficients to be fixed by the initial (and/or
boundary ) conditions.

  * •

calculate the Wronskian of the two independent solutions $y_{1}(x),y_{2}(x)$,

$W(x)=y_{1}(x)y_{2}^{\prime}(x)-y_{1}^{\prime}(x)y_{2}(x).$ (153)

  * •

construct the particular solution of the inhomgeneous equation.

$\begin{split}y_{p}(x)=&y_{1}(x)\int^{x}\frac{y_{2}(z)(-f(z))}{a(z)W(z)}dz+y_{2}(x)\int^{x}\frac{y_{1}(z)f(z)}{a(z)W(z)}dz\\\
=&\int^{x}\frac{y_{1}(z)y_{2}(x)-y_{1}(x)y_{2}(z)}{W(z)}\frac{f(z)}{a(z)}dz.\\\
\end{split}$ (154)

  * •

finally, by adding together the general solution of the homogeneous equation
$y_{h}(x)$ and the particular solution of the inhomogeneous equation
$y_{p}(x)$, obtain the general solution for the inhomogeneous equation

$y(x)=C_{1}y_{1}(x)+C_{2}y_{2}(x)+\int^{x}\frac{y_{1}(z)y_{2}(x)-y_{1}(x)y_{2}(z)}{W(z)}\frac{f(z)}{a(z)}dz.$
(155)

##  Appendix C Integrals for the products of two Bessel functions

There are a few general identities expressing the product of two Bessel
functions in terms of an integral of one Bessel function from Dixion and
Farrar’s paper in 1933 Dixon1933 . For our problem here, only integer orders
of Bessel functions are involved. The more well-known Graf’s formula
Watson1995 ; Abramowitz1965  serves our purpose:

$e^{in\Psi}J_{n}(\omega)=\sum_{m=-\infty}^{\infty}J_{n+m}(Z)J_{m}(Z^{\prime})e^{im\phi}$
(156)

where

$\omega=\sqrt{Z^{2}+Z^{\prime 2}-2ZZ^{\prime}\cos{\phi}}\,.$ (157)

The angle $\Psi$ is defined by

$\omega\cos{\Psi}=Z-Z^{\prime}\cos{\phi},\qquad\omega\sin{\Psi}=Z^{\prime}\sin{\phi}.$
(158)

Integrating both sides of the Graf’s formula using
$\int_{-\pi}^{\pi}\frac{d\phi}{2\pi}e^{-im^{\prime}\phi}$, one obtains

$\int_{-\pi}^{\pi}\frac{d\phi}{2\pi}e^{-im^{\prime}\phi}e^{in\Psi}J_{n}(\omega)=J_{n+m^{\prime}}(Z)J_{m^{\prime}}(Z^{\prime}).$
(159)

This can be rewritten in the form which is more convenient for the
applications

$J_{n+m}(az)J_{m}(bz)=\int_{-\pi}^{\pi}\frac{d\phi}{2\pi}e^{-im\phi}e^{in\Psi}J_{n}(wz)$
(160)

In this case, $w=\sqrt{b^{2}+a^{2}-2ab\cos{\phi}}$ is symmetric with respect
to $a,b$. and

$e^{i\Psi}=\frac{(a-b\cos{\phi})+i(b\sin{\phi})}{\sqrt{b^{2}+a^{2}-2ab\cos{\phi}}}.$
(161)

Here are a few examples that were used in the main body of the paper:

  * $\underline{J_{0}(Z)J_{0}(Z^{\prime})}$:

In this case $m=0,n=0$, thus we obtain

$J_{0}(az)J_{0}(bz)=\int_{-\pi}^{\pi}\frac{d\phi}{2\pi}J_{0}(wz).$ (162)

  * $\underline{J_{0}(Z)J_{1}(Z^{\prime})}$:

Here we can choose between two combinations of $n$ and $m$. For $n=1,m=-1$, we
obtain

$J_{0}(az)J_{1}(bz)=-\int_{-\pi}^{\pi}\frac{d\phi}{2\pi}e^{i\phi}e^{i\Psi}J_{1}(wz)\,.$
(163)

For $m=0,n=1$,

$J_{1}(az)J_{0}(bz)=\int_{-\pi}^{\pi}\frac{d\phi}{2\pi}e^{i\Psi}J_{1}(wz)\,.$
(164)

Therefore we have two integral representation for this product

$J_{0}(az)J_{1}(bz)=-\int_{-\pi}^{\pi}\frac{d\phi}{2\pi}e^{i\phi}e^{i\Psi}J_{1}(wz)=\int_{-\pi}^{\pi}\frac{d\phi}{2\pi}e^{i\Psi^{\prime}}J_{1}(wz)$
(165)

Here $\Psi^{\prime}$ is related to $\Psi$ by the exchange $a\leftrightarrow
b$,

$e^{i\Psi^{\prime}}=\frac{(b-a\cos{\phi})+i(a\sin{\phi})}{\sqrt{b^{2}+a^{2}-2ab\cos{\phi}}}.$
(166)

It can be easily checked that $\phi+\pi+\Psi=\Psi^{\prime}$. Thus the two
integral representations are equivalent. Additionally, since $w$ is an even
function of $\phi$, only $\phi-$even part of the exponentials contributes:

$\begin{split}J_{0}(az)J_{1}(bz)=&2\int_{0}^{\pi}\frac{d\phi}{2\pi}\cos{\Psi^{\prime}}J_{1}(wz)\\\
=&2\int_{|b-a|}^{|b+a|}dw\frac{b^{2}-a^{2}+w^{2}}{b\sqrt{((a+b)^{2}-w^{2})(w^{2}-(a-b)^{2})}}J_{1}(wz).\end{split}$
(167)

In the second equality, we have changed the integration variable from $\phi$
to $w$ using

$\cos\phi=\frac{a^{2}+b^{2}-w^{2}}{2ab},$ (168)
$d\phi=\frac{wdw}{ab\sin{\phi}}=\frac{2wdw}{\sqrt{4a^{2}b^{2}-(a^{2}+b^{2}-w^{2})^{2}}}.$
(169)
$\cos\Psi^{\prime}=\frac{b-a\cos{\phi}}{w}=\frac{b^{2}-a^{2}+w^{2}}{2bw}.$
(170)

  * $\underline{J_{1}(Z)J_{1}(Z^{\prime})}$:

Setting $n=0,m=1$, we get

$J_{1}(az)J_{1}(bz)=\int_{-\pi}^{\pi}\frac{d\phi}{2\pi}e^{-i\phi}J_{0}(wz).$
(171)

Alternatively for $m=-1,n=0$,

$J_{1}(az)J_{1}(bz)=\int_{-\pi}^{\pi}\frac{d\phi}{2\pi}e^{i\phi}J_{0}(wz).$
(172)

These two expressions are equivalent. Since $w$ is even function of $\phi$,
only the even part of $e^{i\phi}$ contributes to the integral:

$J_{1}(az)J_{1}(bz)=2\int_{0}^{\pi}\frac{d\phi}{2\pi}\cos{\phi}J_{0}(wz).$
(173)

On the other hand, for $n=2,m=-1$,

$J_{1}(az)J_{1}(bz)=-\int_{-\pi}^{\pi}\frac{d\phi}{2\pi}e^{i\phi}e^{i2\Psi}J_{2}(wz).$
(174)

The case of $n=-2,m=1$ will yield the same expression. We can further simplify
to get

$J_{1}(az)J_{1}(bz)=-2\int_{0}^{\pi}\frac{d\phi}{2\pi}\cos{(\phi+2\Psi)}J_{2}(wz).$
(175)

  * $\underline{J_{1}(Z)J_{2}(Z^{\prime})}$:

There are two ways to express $J_{1}(Z)J_{2}(Z^{\prime})$ in terms of an
integral of $J_{1}(Z^{\prime\prime})$. For $n=1,m=-2$ or $n=-1,m=2$, one
obtains

$J_{1}(az)J_{2}(bz)=-\int_{-\pi}^{\pi}\frac{d\phi}{2\pi}e^{2i\phi}e^{i\Psi}J_{1}(wx).$
(176)

For $n=1,m=1$ or $n=-1,m=-1$,

$J_{2}(az)J_{1}(bz)=\int_{-\pi}^{\pi}\frac{d\phi}{2\pi}e^{-i\phi}e^{i\Psi}J_{1}(wx).$
(177)

Again, it seems like we have two expressions

$J_{1}(az)J_{2}(bz)=-\int_{-\pi}^{\pi}\frac{d\phi}{2\pi}e^{2i\phi}e^{i\Psi}J_{1}(wx)=\int_{-\pi}^{\pi}\frac{d\phi}{2\pi}e^{-i\phi}e^{i\Psi^{\prime}}J_{1}(wx),$
(178)

but they are equivalent due to
$\cos{(\phi-\Psi^{\prime})}=-\cos{(2\phi+\Psi)}$.

##  Appendix D Solutions in the non-Coulomb subgauge

In this appendix, we present the solutions of the Yang-Mills equation at
order-$g$ and order-$g^{3}$ in the non-Coulomb sub-gauge that is defined by
$U(\mathbf{x})$. Initial conditions in this gauge has been discussed in Sec.
3.2.1.

Order-$g$ solutions

$\begin{split}&\zeta^{(1)}(\tau,\mathbf{k})=b_{\eta}(\mathbf{k})\frac{J_{1}(k_{\perp}\tau)}{k_{\perp}\tau},\\\
&\zeta^{(1)}_{i}(\tau,\mathbf{k})=\frac{-i\epsilon^{il}\mathbf{k}_{l}}{k_{\perp}^{2}}b_{\perp}(\mathbf{k})J_{0}(k_{\perp}\tau)-i\mathbf{k}_{i}\Lambda(\mathbf{k})\\\
\end{split}$ (179)

with the initial conditions

$\begin{split}&b_{\eta}(\mathbf{k})=2\zeta^{(1)}(\tau=0,\mathbf{k}),\\\
&b_{\perp}(\mathbf{k})=-i\epsilon^{ij}\mathbf{k}_{i}\zeta_{j}^{(1)}(\tau=0,\mathbf{k}),\\\
&\Lambda(\mathbf{k})=\frac{ik_{i}}{k^{2}}\zeta_{i}^{(1)}(\tau=0,\mathbf{k}).\\\
\end{split}$ (180)

Here $b_{\eta}(\mathbf{k})$ and $b_{\perp}(\mathbf{k})$ are exactly the same
as those in Eq. (53). In the non-Coulomb subgauge, $\partial_{i}\zeta_{i}\neq
0$, so $\Lambda(\mathbf{k})$ contributes to the initial condition.

The order-$g^{3}$ solutions are

$\begin{split}\zeta^{(3)}(\tau,\mathbf{k})=&2\zeta^{(3)}(\tau=0,\mathbf{k})\frac{J_{1}(k_{\perp}\tau)}{k_{\perp}\tau}+i\int\frac{d^{2}\mathbf{p}}{(2\pi)}\Bigg{(}\Big{[}\Lambda(\mathbf{k}-\mathbf{p}),b_{\eta}(\mathbf{p})\Big{]}\left(\frac{J_{1}(p_{\perp}\tau)}{p_{\perp}\tau}-\frac{J_{1}(k_{\perp}\tau)}{k_{\perp}\tau}\right)\\\
&-\frac{\mathbf{k}\times\mathbf{p}}{p_{\perp}^{2}|\mathbf{k}-\mathbf{p}|^{2}}\Big{[}b_{\perp}(\mathbf{k}-\mathbf{p}),b_{\eta}(\mathbf{p})\Big{]}\int_{-\pi}^{\pi}\frac{d\phi}{2\pi}\left(-1+\frac{2\mathbf{k}\cdot\mathbf{p}}{k_{\perp}^{2}-w_{\perp}^{2}}\right)\left(\frac{J_{1}(w_{\perp}\tau)}{w_{\perp}\tau}-\frac{J_{1}(k_{\perp}\tau)}{k_{\perp}\tau}\right)\Bigg{)},\end{split}$
(181)
$\begin{split}\zeta^{(3)}_{\perp}(\tau,\mathbf{k})=&\zeta^{(3)}_{\perp}(\tau=0,\mathbf{k})J_{0}(k_{\perp}\tau)-\frac{i}{k_{\perp}}\int\frac{d^{2}\mathbf{p}}{(2\pi)^{2}}\int_{-\pi}^{\pi}\frac{d\phi}{2\pi}\frac{1}{k_{\perp}^{2}-w_{\perp}^{2}}\\\
&\times\frac{(\mathbf{k}\times\mathbf{p})}{2p^{2}_{\perp}|\mathbf{k}-\mathbf{p}|^{2}}\Big{(}(p_{\perp}^{2}+|\mathbf{k}-\mathbf{p}|^{2}-w_{\perp}^{2})\Big{[}b_{\eta}(\mathbf{p}),b_{\eta}(\mathbf{k}-\mathbf{p})\Big{]}\\\
&+2(\mathbf{p}\cdot\mathbf{k}-p_{\perp}^{2}-k_{\perp}^{2})\Big{[}b_{\perp}(\mathbf{p}),b_{\perp}(\mathbf{k}-\mathbf{p})\Big{]}\Big{)}\left(J_{0}(w_{\perp}\tau)-J_{0}(k_{\perp}\tau)\right)\\\
&+\frac{(\mathbf{k}\cdot\mathbf{p})}{p_{\perp}^{2}}\Big{[}b_{\perp}(\mathbf{p}),\Lambda(\mathbf{k}-\mathbf{p})\Big{]}(J_{0}(p_{\perp}\tau)-J_{0}(k_{\perp}\tau))\\\
&+\frac{1}{2}(\mathbf{p}\times\mathbf{k})\Big{[}\Lambda(\mathbf{p}),\Lambda(\mathbf{k}-\mathbf{p})\Big{]}(1-J_{0}(k_{\perp}\tau)),\\\
\end{split}$ (182)
$\begin{split}\Lambda^{(3)}(\tau,\mathbf{k})=&\Lambda^{(3)}(\tau=0,\mathbf{k})-\frac{i}{k_{\perp}^{2}}\int\frac{d^{2}\mathbf{p}}{(2\pi)^{2}}\int_{-\pi}^{\pi}\frac{d\phi}{2\pi}\frac{\mathbf{k}\cdot(\mathbf{k}-2\mathbf{p})}{4w^{2}_{\perp}p_{\perp}^{2}|\mathbf{k}-\mathbf{p}|^{2}}\\\
&\Big{(}2\mathbf{p}\cdot(\mathbf{k}-\mathbf{p})\left[b_{\perp}(\mathbf{p}),b_{\perp}(\mathbf{k}-\mathbf{p})\right]-(k_{\perp}^{2}-w_{\perp}^{2}-2\mathbf{p}(\mathbf{k}-\mathbf{p}))\left[b_{\eta}(\mathbf{p}),b_{\eta}(\mathbf{k}-\mathbf{p})\right]\Big{)}\\\
&\times(1-J_{0}(w_{\perp}\tau))+\frac{\mathbf{p}\times\mathbf{k}}{|\mathbf{k}-\mathbf{p}|^{2}}\left[\Lambda(\mathbf{p}),b_{\perp}(\mathbf{k}-\mathbf{p})\right](1-J_{0}(|\mathbf{k}-\mathbf{p}|\tau)).\\\
\end{split}$ (183)

Comparing with the solutions in the initial time Coulomb sub gauge given by
Eqs. (70), (83) and (95), the terms containing $\Lambda(\mathbf{p})$ in Eqs.
(181), (182), (183), which represent final state interactions in the non-
Coulomb sub gague, are shifted to become initial state effects in Eqs. (70),
(83) and (95)

## References

  * (1) C. Gale, S. Jeon and B. Schenke, _Hydrodynamic Modeling of Heavy-Ion Collisions_ , _Int. J. Mod. Phys. A_ 28 (2013) 1340011 [1301.5893].
  * (2) U. Heinz and R. Snellings, _Collective flow and viscosity in relativistic heavy-ion collisions_ , _Ann. Rev. Nucl. Part. Sci._ 63 (2013) 123 [1301.2826].
  * (3) B. Schenke, C. Shen and P. Tribedy, _Hybrid Color Glass Condensate and hydrodynamic description of the Relativistic Heavy Ion Collider small system scan_ , _Phys. Lett. B_ 803 (2020) 135322 [1908.06212].
  * (4) R. Baier, A.H. Mueller, D. Schiff and D.T. Son, _’Bottom up’ thermalization in heavy ion collisions_ , _Phys. Lett. B_ 502 (2001) 51 [hep-ph/0009237].
  * (5) A. Kurkela and Y. Zhu, _Isotropization and hydrodynamization in weakly coupled heavy-ion collisions_ , _Phys. Rev. Lett._ 115 (2015) 182301 [1506.06647].
  * (6) A. Kurkela, A. Mazeliauskas, J.-F. Paquet, S. Schlichting and D. Teaney, _Matching the Nonequilibrium Initial Stage of Heavy Ion Collisions to Hydrodynamics with QCD Kinetic Theory_ , _Phys. Rev. Lett._ 122 (2019) 122302 [1805.01604].
  * (7) B. Wu and Y.V. Kovchegov, _Time-dependent observables in heavy ion collisions. Part I. Setting up the formalism_ , _JHEP_ 03 (2018) 158 [1709.02866].
  * (8) Y.V. Kovchegov and B. Wu, _Time-dependent observables in heavy ion collisions. Part II. In search of pressure isotropization in the $\varphi^{4}$ theory_, _JHEP_ 03 (2018) 157 [1709.02868].
  * (9) A. Krasnitz and R. Venugopalan, _The Initial energy density of gluons produced in very high-energy nuclear collisions_ , _Phys. Rev. Lett._ 84 (2000) 4309 [hep-ph/9909203].
  * (10) A. Krasnitz and R. Venugopalan, _The Initial gluon multiplicity in heavy ion collisions_ , _Phys. Rev. Lett._ 86 (2001) 1717 [hep-ph/0007108].
  * (11) B. Schenke, P. Tribedy and R. Venugopalan, _Fluctuating Glasma initial conditions and flow in heavy ion collisions_ , _Phys. Rev. Lett._ 108 (2012) 252301 [1202.6646].
  * (12) C. Gale, S. Jeon, B. Schenke, P. Tribedy and R. Venugopalan, _Event-by-event anisotropic flow in heavy-ion collisions from combined Yang-Mills and viscous fluid dynamics_ , _Phys. Rev. Lett._ 110 (2013) 012302 [1209.6330].
  * (13) K. Dusling, W. Li and B. Schenke, _Novel collective phenomena in high-energy proton–proton and proton–nucleus collisions_ , _Int. J. Mod. Phys. E_ 25 (2016) 1630002 [1509.07939].
  * (14) J.L. Nagle and W.A. Zajc, _Small System Collectivity in Relativistic Hadronic and Nuclear Collisions_ , _Ann. Rev. Nucl. Part. Sci._ 68 (2018) 211 [1801.03477].
  * (15) Y.V. Kovchegov and A.H. Mueller, _Gluon production in current nucleus and nucleon - nucleus collisions in a quasiclassical approximation_ , _Nucl. Phys. B_ 529 (1998) 451 [hep-ph/9802440].
  * (16) B.Z. Kopeliovich, A.V. Tarasov and A. Schafer, _Bremsstrahlung of a quark propagating through a nucleus_ , _Phys. Rev. C_ 59 (1999) 1609 [hep-ph/9808378].
  * (17) A. Kovner and U.A. Wiedemann, _Eikonal evolution and gluon radiation_ , _Phys. Rev. D_ 64 (2001) 114002 [hep-ph/0106240].
  * (18) A. Dumitru and L.D. McLerran, _How protons shatter colored glass_ , _Nucl. Phys. A_ 700 (2002) 492 [hep-ph/0105268].
  * (19) J.P. Blaizot, F. Gelis and R. Venugopalan, _High energy pA collisions in the color glass condensate approach I: Gluon production and the Cronin effect_ , _Nucl. Phys._ A743 (2004) 13 [hep-ph/0402256].
  * (20) I. Balitsky, _Scattering of shock waves in QCD_ , _Phys. Rev. D_ 70 (2004) 114030 [hep-ph/0409314].
  * (21) G.A. Chirilli, Y.V. Kovchegov and D.E. Wertepny, _Classical Gluon Production Amplitude for Nucleus-Nucleus Collisions: First Saturation Correction in the Projectile_ , _JHEP_ 03 (2015) 015 [1501.03106].
  * (22) L. McLerran and V. Skokov, _Odd Azimuthal Anisotropy of the Glasma for pA Scattering_ , _Nucl. Phys. A_ 959 (2017) 83 [1611.09870].
  * (23) Y.V. Kovchegov and V.V. Skokov, _How classical gluon fields generate odd azimuthal harmonics for the two-gluon correlation function in high-energy collisions_ , _Phys. Rev. D_ 97 (2018) 094021 [1802.08166].
  * (24) I. Balitsky, _Operator expansion for high-energy scattering_ , _Nucl. Phys._ B463 (1996) 99 [hep-ph/9509348].
  * (25) I. Balitsky, _Factorization and high-energy effective action_ , _Phys. Rev._ D60 (1999) 014020 [hep-ph/9812311].
  * (26) I. Balitsky, _Effective field theory for the small-x evolution_ , _Phys. Lett._ B518 (2001) 235 [hep-ph/0105334].
  * (27) J. Jalilian-Marian, A. Kovner, A. Leonidov and H. Weigert, _The BFKL equation from the Wilson renormalization group_ , _Nucl. Phys._ B504 (1997) 415 [hep-ph/9701284].
  * (28) J. Jalilian-Marian, A. Kovner, A. Leonidov and H. Weigert, _The Wilson renormalization group for low x physics: Towards the high density regime_ , _Phys. Rev._ D59 (1999) 014014 [hep-ph/9706377].
  * (29) J. Jalilian-Marian, A. Kovner and H. Weigert, _The Wilson renormalization group for low x physics: Gluon evolution at finite parton density_ , _Phys. Rev._ D59 (1999) 014015 [hep-ph/9709432].
  * (30) E. Iancu, A. Leonidov and L.D. McLerran, _The renormalization group equation for the color glass condensate_ , _Phys. Lett._ B510 (2001) 133 [hep-ph/0102009].
  * (31) E. Iancu, A. Leonidov and L.D. McLerran, _Nonlinear gluon evolution in the color glass condensate. i_ , _Nucl. Phys._ A692 (2001) 583 [hep-ph/0011241].
  * (32) E. Ferreiro, E. Iancu, A. Leonidov and L. McLerran, _Nonlinear gluon evolution in the color glass condensate. ii_ , _Nucl. Phys._ A703 (2002) 489 [hep-ph/0109115].
  * (33) H. Weigert, _Unitarity at small Bjorken x_ , _Nucl. Phys._ A703 (2002) 823 [hep-ph/0004044].
  * (34) L.D. McLerran and R. Venugopalan, _Computing quark and gluon distribution functions for very large nuclei_ , _Phys. Rev. D_ 49 (1994) 2233 [hep-ph/9309289].
  * (35) L.D. McLerran and R. Venugopalan, _Gluon distribution functions for very large nuclei at small transverse momentum_ , _Phys. Rev. D_ 49 (1994) 3352 [hep-ph/9311205].
  * (36) T. Lappi, _Wilson line correlator in the MV model: Relating the glasma to deep inelastic scattering_ , _Eur. Phys. J. C_ 55 (2008) 285 [0711.3039].
  * (37) Y.V. Kovchegov, _Quantum structure of the nonAbelian Weizsacker-Williams field for a very large nucleus_ , _Phys. Rev. D_ 55 (1997) 5445 [hep-ph/9701229].
  * (38) S. Schlichting and V. Skokov, _Saturation corrections to dilute-dense particle production and azimuthal correlations in the Color Glass Condensate_ , _Phys. Lett. B_ 806 (2020) 135511 [1910.12496].
  * (39) A. Kovner and M. Lublinsky, _Angular and long range rapidity correlations in particle production at high energy_ , _Int. J. Mod. Phys. E_ 22 (2013) 1330001 [1211.1928].
  * (40) Y.V. Kovchegov and D.E. Wertepny, _Long-Range Rapidity Correlations in Heavy-Light Ion Collisions_ , _Nucl. Phys. A_ 906 (2013) 50 [1212.1195].
  * (41) I. Balitsky, _High-energy effective action from scattering of QCD shock waves_ , _Phys. Rev. D_ 72 (2005) 074027 [hep-ph/0507237].
  * (42) Y. Hatta, E. Iancu, L. McLerran, A. Stasto and D. Triantafyllopoulos, _Effective Hamiltonian for QCD evolution at high energy_ , _Nucl. Phys. A_ 764 (2006) 423 [hep-ph/0504182].
  * (43) A. Kovner, E. Levin, M. Li and M. Lublinsky, _Reggeon Field Theory and Self Duality: Making Ends Meet_ , _JHEP_ 10 (2020) 185 [2007.12132].
  * (44) A. Krasnitz, Y. Nara and R. Venugopalan, _Gluon production in the color glass condensate model of collisions of ultrarelativistic finite nuclei_ , _Nucl. Phys. A_ 717 (2003) 268 [hep-ph/0209269].
  * (45) T. Lappi, _Production of gluons in the classical field model for heavy ion collisions_ , _Phys. Rev. C_ 67 (2003) 054903 [hep-ph/0303076].
  * (46) B. Schenke, S. Schlichting and R. Venugopalan, _Azimuthal anisotropies in p $+$Pb collisions from classical Yang–Mills dynamics_, _Phys. Lett. B_ 747 (2015) 76 [1502.01331].
  * (47) L.D. McLerran, _The Color glass condensate and small x physics: Four lectures_ , _Lect. Notes Phys._ 583 (2002) 291 [hep-ph/0104285].
  * (48) E. Iancu and R. Venugopalan, _The Color glass condensate and high-energy scattering in QCD_ , in _Quark-gluon plasma 4_ , R.C. Hwa and X.-N. Wang, eds., pp. 249–3363 (2003), DOI [hep-ph/0303204].
  * (49) Y.V. Kovchegov and E. Levin, _Quantum chromodynamics at high energy_ , vol. 33, Cambridge University Press (8, 2012), 10.1017/CBO9781139022187.
  * (50) A. Kovner, L.D. McLerran and H. Weigert, _Gluon production from nonAbelian Weizsacker-Williams fields in nucleus-nucleus collisions_ , _Phys. Rev. D_ 52 (1995) 6231 [hep-ph/9502289].
  * (51) M. Gyulassy and L.D. McLerran, _Yang-Mills radiation in ultrarelativistic nuclear collisions_ , _Phys. Rev. C_ 56 (1997) 2219 [nucl-th/9704034].
  * (52) J.-P. Blaizot, T. Lappi and Y. Mehtar-Tani, _On the gluon spectrum in the glasma_ , _Nucl. Phys. A_ 846 (2010) 63 [1005.0955].
  * (53) J. Berges, K. Boguslavski, S. Schlichting and R. Venugopalan, _Universal attractor in a highly occupied non-Abelian plasma_ , _Phys. Rev. D_ 89 (2014) 114007 [1311.3005].
  * (54) M. Li and V. Skokov, “First Saturation Correction in High Energy Proton-Nucleus Collisions: II. Single Inclusive Soft Gluon Production.” In preparation.
  * (55) R. Fries, J. Kapusta and Y. Li, _Near-fields and initial energy density in the color glass condensate model_ , nucl-th/0604054.
  * (56) G. Chen, R.J. Fries, J.I. Kapusta and Y. Li, _Early Time Dynamics of Gluon Fields in High Energy Nuclear Collisions_ , _Phys. Rev. C_ 92 (2015) 064912 [1507.03524].
  * (57) A.L. Dixon and W.L. Ferrar, _Integrals for the Product of Two Bessel Functions (II)_ , _The Quarterly Journal of Mathematics_ os-4 (1933) 297.
  * (58) G.N. Watson, _A Treatise on the Theory of Bessel Functions_ , 2nd edition, Cambridge University Press (1995).
  * (59) M. Abramowitz and I. Stegun, _Handbook of Mathematical Functions: with Formulas, Graphs, and Mathematical Tables_ , 0009-Revised edition, Dover Publications (1965).