11institutetext: University of Tübingen 22institutetext: University of Lucerne
22email<EMAIL_ADDRESS>
# Segmentation-guided MRI reconstruction for meaningfully diverse
reconstructions

Jan Nikolas Morshuis 11 Matthias Hein 11 Christian F. Baumgartner 1122

###### Abstract

Inverse problems, such as accelerated MRI reconstruction, are ill-posed and an
infinite amount of possible and plausible solutions exist. This may not only
lead to uncertainty in the reconstructed image but also in downstream tasks
such as semantic segmentation. This uncertainty, however, is mostly not
analyzed in the literature, even though probabilistic reconstruction models
are commonly used. These models can be prone to ignore plausible but unlikely
solutions like rare pathologies. Building on MRI reconstruction approaches
based on diffusion models, we add guidance to the diffusion process during
inference, generating two meaningfully diverse reconstructions corresponding
to an upper and lower bound segmentation. The reconstruction uncertainty can
then be quantified by the difference between these bounds, which we coin the
’uncertainty boundary’. We analyzed the behavior of the upper and lower bound
segmentations for a wide range of acceleration factors and found the
uncertainty boundary to be both more reliable and more accurate compared to
repeated sampling. Code is available at
https://github.com/NikolasMorshuis/SGR.

##  1 Introduction

Accelerated MRI reconstruction has substantially improved in recent years
through the use of machine learning methods. Currently, diffusion models (DMs)
[13, 19] achieve state-of-the-art MRI reconstruction scores [4, 14, 5]. An
often neglected area in the literature is the inherent uncertainty of the
reconstruction process. Accelerated MRI reconstruction is an ill-posed problem
and infinite plausible reconstructions exist for an undersampled MRI image.
However, in the majority of recent work only a single point estimate (a single
reconstruction) of the solution space is provided. This common practice is
problematic in safety-critical applications like medical imaging, as it does
not accurately reflect the inherent uncertainty and no estimation of the
confidence of the solution can be provided.

Figure 1: Motivation. Repeated reconstruction tends to create MRI
reconstructions that are conceptually similar and not meaningfully diverse.
Similar results are also observed in downstream tasks like segmentation in our
experiments. Our method in contrast only reconstructs two images, an upper
bound ($x_{0}^{\uparrow}$) and a lower bound ($x_{0}^{\downarrow}$)
reconstruction, corresponding to an upper ($V_{seg}^{\uparrow}$) and lower
bound ($V_{seg}^{\downarrow}$) on the segmentation volume, giving an intuitive
understanding of the segmentation uncertainty.

A common approach to quantify the reconstruction uncertainty is to reconstruct
a large number of samples and analyze their pixel-wise variance. This kind of
uncertainty prediction has been investigated with probabilistic networks [12,
3], Monte Carlo dropout [16] ensembling [15] or DMs [5, 14, 24].

Recently however, [6] has shown that Repeated Reconstruction (RR) in natural
image reconstruction with DMs leads to similar images which display a lack of
meaningful diversity. Even though the diffusion model is trained to model the
posterior distribution accurately, the heavy tail of the true posterior
distribution is often largely ignored in practice and plausible yet unlikely
solutions with different semantic meaning are not found. This is a problem for
medical imaging, where certain pathologies might be rare and unlikely.
However, it is important to recover these cases precisely and include them as
possible solutions when reconstructing multiple samples for uncertainty
prediction.

The definition of meaningful diversity can depend on the downstream task and
might be different for e.g. pathology detection (a pathology exists or not) or
semantic segmentation (the segmented tissue might be smaller or larger). In
this work we focus on diversity of semantic segmentation, as it is an
important task in clinical practice and the results are more easily
interpretable for non-experts.

In this paper, instead of random sampling, we propose to reconstruct two MRI-
images per segmentation class corresponding to the upper and lower bound
segmentation volume. To this end, we introduce a novel method called
Segmentation Guided Reconstruction (SGR). SGR builds on DM-based MRI
reconstruction. However, rather than standard sampling, we guide the
reconstruction process with the gradient of segmentation losses that either
penalize small segmentation volumes or large segmentation volumes. The
resulting lower and upper bound reconstructions and segmentations are an
intuitive way of understanding the segmentation uncertainty which arises due
to the ill-posed nature of the reconstruction problem. See Fig. 1 for a visual
explanation of the motivation. In summary, our contributions are:

  * •

We introduce a novel method Segmentation-Guided Reconstruction (SGR) to guide
the reconstruction process of undersampled MRI data, leading to meaningfully
diverse reconstructions.

  * •

We show that the diverse reconstructions lead to diverse segmentations that
correspond to an upper and lower bound on the segmentation. We introduce the
concept of ’uncertainty boundary’ for the volume between both bounds.

  * •

We show that our method captures the inherent uncertainty more faithfully than
the commonly used Repeated Reconstruction (RR), potentially leading to safer
MRI reconstruction methods

##  2 Related Work

Diffusion based inverse problem solvers have seen large success in MRI
reconstruction [14, 5, 24, 18], achieving state-of-the-art reconstruction
scores [14, 5, 4]. In contrast to supervised learning-based methods [20, 25],
trained DMs can be used for several acceleration factors. Recent
reconstruction methods [4] allow for more than $80\times$ faster inference
time, reducing the computational burden for RR. With RR it is possible to
reconstruct several distinct images, allowing to estimate the reconstruction
uncertainty by analyzing the variance of the voxel values. However, the
distinct images do not necessarily display meaningful diversity [6].

While the reconstruction and segmentation problem are mostly tackled
separately, some work exists aiming to improve segmentation results by
combining reconstruction and segmentation in an end-to-end fashion [2, 21, 1,
23, 22]. Wu et al. [23] show that fine-tuning a pre-trained reconstruction
network with a task-specific network head leads to improved segmentation
results. Similarly, the winning method of the K2S challenge [22] also first
trains a reconstruction network before fine-tuning the network for
segmentation. These approaches, however, lead to degraded reconstructions.
Moreover, all of these methods lack an analysis of the effects of the inherent
reconstruction uncertainty on segmentation.

The literature on reconstruction uncertainty propagation to downstream tasks
is relatively sparse and most methods utilize repeated sampling of
probabilistic networks that might not offer meaningful diversity [6]. Feiner
et al. [11] propose to utilize Monte Carlo sampling to estimate the
reconstruction uncertainty that can then be propagated to classification
tasks. Cheung et al. [3] use repeated sampling of CT reconstructions to obtain
a distribution of several downstream metrics, which they use to calculate
statistically valid prediction sets using conformal prediction. Fischer et al.
[12] are the first to consider the propagation of the aleatoric uncertainty
due to the reconstruction process in the context of segmentation. By repeated
sampling of a probabilistic reconstruction network and subsequent
segmentation, they were able to quantify the variance of the segmentation at
pixel-level. However, data-consistency was not enforced and no upper or lower
bound of the segmentation has been estimated.

##  3 Method

###  3.1 Background diffusion models and MRI-reconstruction

In the inverse problem of MRI reconstruction, we aim to reconstruct a
plausible image $\bm{x}_{0}$ that is coherent with the measurement data
$\bm{y}$ when passing through a forward operator (i.e. Fourier Transform +
Masking) $\bm{A}$, such that the following equation holds:

$\bm{y}=\bm{A}\bm{x}_{0}^{(i)}\ \text{.}$ (1)

The problem is typically ill-posed and infinite $\bm{x}_{0}^{(i)}$ exist that
are possible solutions for Eq. 1. In recent years DMs have shown to perform
exceptionally well on the MRI-reconstruction tasks, being able to sample a set
of distinct solutions $\bm{x}_{0}^{(i)}$ from the posterior.

DMs generate samples from the data distribution $p_{data}(\bm{x}_{0})$ by
reversing a diffusion process from timestep $t=0$ to $t=T$ described by a
conditional density

$q(\bm{x}_{t}|\bm{x}_{0})=\mathcal{N}(\bm{x}_{t}|\sqrt{\bar{\alpha_{t}}}\bm{x}_{0},(1-\bar{\alpha_{t}})\bm{I})$
(2)

with a pre-defined defined $\bar{\alpha_{t}}$ schedule.

A DM can then be trained by predicting the noise through epsilon matching:

$\min\limits_{\theta}\mathbb{E}_{\bm{x}_{t}\sim
q(\bm{x}_{t}|\bm{x}_{0}),\bm{x}_{0}\sim
p_{data}(\bm{x}_{0}),\bm{\epsilon}\sim\mathcal{N}(0,\bm{I})}(||\bm{\epsilon}_{\theta}^{(t)}(\bm{x}_{t})-\bm{\epsilon}||_{2}^{2})\
\text{.}$ (3)

The learned $\bm{\epsilon}$-matching function can then be used to predict the
clean image $\hat{\bm{x}}_{0|t}$ at every step $t$ of the diffusion process
using Tweedie denoising [10]:

$\hat{\bm{x}}_{0|t}=(\bm{x}_{t}-\sqrt{1-\bar{\alpha_{t}}}\bm{\epsilon}_{\theta}^{(t)}(\bm{x}_{t}))/\sqrt{\bar{\alpha_{t}}}\text{.}$
(4)

Recently, Chung et al. [4] introduced a DM based reconstruction method that
utilizes conjugate gradient (CG) optimization: During the CG-optimization
process, starting from $\hat{\bm{x}}_{0|t}$ we search a
$\hat{\bm{x}}_{0|t}^{\prime}$ that minimizes the distance
$\bm{A}^{*}\bm{A}\hat{\bm{x}}_{0|t}^{\prime}-\bm{A}^{*}\bm{y}$ with
$\bm{A}^{*}$ being the backward operator mapping from the measurement space
back to the image space:

$\hat{\bm{x}}_{0|t}^{\prime}=CG(\bm{A}^{*}\bm{A},\bm{A}^{*}\bm{y},\hat{\bm{x}}_{0|t})\text{.}$
(5)

The update rule using the DDIM method [17] is then given as follows:

$\bm{x}_{t-1}=\sqrt{\bar{\alpha}_{t-1}}\hat{\bm{x}}_{0|t}^{\prime}+\sqrt{1-\bar{\alpha}_{t-1}-\eta^{2}\tilde{\beta}_{t}^{2}}\bm{\epsilon}_{\theta}^{(t)}(x_{t})+\eta\tilde{\beta}_{t}\bm{\epsilon}\text{.}$
(6)

We modify the algorithm from [4] slightly to achieve full data-consistency.
For the last diffusion step $t=1\rightarrow t=0$, we replace Eq. 5 and make
sure that the reconstructed image $\bm{x}_{0}$ is fully data-consistent using
$\bm{y}$ and Eq. 1.

###  3.2 Segmentation Guidance for diverse sampling

Figure 2: Method explanation. For each inverse diffusion step $t$, we
calculate $\nabla_{x_{t}}\mathcal{L}^{\uparrow/\downarrow}$ and include the
gradient in the calculation of $x_{t-1}$, in order to increase (decrease) the
segmentation volume and to get an upper- and lower-bound segmentation
($S^{\uparrow}$ and $S^{\downarrow}$).

In order to generate a set of meaningfully diverse solutions
$\bm{x}_{0}^{(i)}$, we reconstruct two images per segmentation class
($\bm{x}_{0}^{\uparrow}$ and $\bm{x}_{0}^{\downarrow}$) with corresponding
upper and a lower bound segmentations ($S^{\uparrow}$ and $S^{\downarrow}$) as
well as segmentation volumes ($V_{seg}^{\uparrow}$ and
$V_{seg}^{\downarrow}$). The volume between $S^{\downarrow}$ and
$S^{\uparrow}$ is the uncertainty boundary $S_{unc}$. We guide the
reconstruction process of the diffusion model using segmentation guidance. To
achieve this, we propose two loss-functions $\mathcal{L}^{\uparrow}$ and
$\mathcal{L}^{\downarrow}$ that are suitable for segmentation-based model
guidance.

Let $p_{\phi}^{c}(\bm{x}_{0,jk})$ be the predicted probability of pixel
$(j,k)$ in the image $\bm{x}_{0}$ belonging to the segmentation class
$c\in\mathcal{C}$, and let $p_{\phi}$ be a segmentation network. We abbreviate
the term as $\bm{p}_{jk}^{(c)}$. The volume of $c$ can be calculated by
summing the number of pixels that belong to $c$ multiplied by voxel-volume
$V_{voxel}$:

$V_{c}(\bm{x}_{0})=\sum_{j,k}(\mathbb{I}_{c}(\bm{p}_{jk}^{(c)}))\cdot
V_{voxel}\text{,}$ (7)

where $\mathbb{I}_{c}$ is an indicator function

$\mathbb{I}_{c}=\begin{cases}1&\text{if
}\bm{p}_{jk}^{(c)}>\bm{p}_{jk}^{(c^{\prime})}\text{ }\forall
c^{\prime}\in\mathcal{C},c^{\prime}\neq c\\\
0&$\text{otherwise}$\end{cases}\text{.}$ (8)

During optimization, we consider only one class at a time, such that we can
simplify the objective to a binary segmentation class. We build our guidance-
loss based on the widely-used Binary-Cross-Entropy loss:

$\mathcal{L}_{BCE}(\bm{y},\bm{p},c)=-\sum_{j,k}\Bigl{(}\bm{y}_{jk}\log(\bm{p}_{jk}^{(c)})+(1-\bm{y}_{jk})\log(1-\bm{p}_{jk}^{(c)})\Bigr{)}\text{.}$
(9)

We can maximize the optimization objective $\mathcal{L}^{\uparrow}$ to find
the upper bound segmentation, as well as $\mathcal{L}^{\downarrow}$ to find
the lower bound segmentation:

$\mathcal{L}^{\uparrow}=\mathcal{L}_{BCE}(\bm{0},\bm{p},c)=-\sum_{j,k}\log(1-\bm{p}_{jk}^{(c)})\text{,}$
(10)

$\mathcal{L}^{\downarrow}=\mathcal{L}_{BCE}(\bm{1},\bm{p},c)=-\sum_{j,k}\log(\bm{p}_{jk}^{(c)})\text{.}$
(11)

The binary cross-entropy loss is unbounded, meaning that pixel probabilities
are further optimized, even if the target class has already been reached.
Following Croce et al. [7] we can avoid this behavior by using a masked cross-
entropy loss, which excludes pixels already segmented as the target class from
the optimization process. The optimization objectives can therefore be written
as follows:

$\mathcal{L}_{m}^{\downarrow}=\mathbb{I}_{c}\mathcal{L}^{\downarrow};\quad\mathcal{L}_{m}^{\uparrow}=(1-\mathbb{I}_{c})\mathcal{L}^{\uparrow}\text{.}$
(12)

Figure 3: Sampling vs. adversarial guidance. Examples of lower-
($x_{0}^{\downarrow}$) and upper-bound ($x_{0}^{\uparrow}$) reconstructions
(16x acc.) and segmentations ($S^{\downarrow}$, $S^{\uparrow}$) using our SGR
method and the standard RR method. ($S^{f}$ is the segmentation of the fully-
sampled image. Green: $S^{f}\&S^{\uparrow}$ or $S^{f}\&S^{\downarrow}$, Blue:
$S^{f}>S^{\downarrow}$ or $S^{f}<S^{\uparrow}$, Red: $S^{f}<S^{\downarrow}$ or
$S^{f}>S^{\uparrow}$)

The derived optimization objective can be used during the reverse diffusion
process. We therefore calculate the gradient of $\mathcal{L}_{m}^{\uparrow}$
(or $\mathcal{L}_{m}^{\downarrow}$, respectively) with respect to $x_{t}$ at
every step $t$, to find the perturbation necessary to increase the
optimization loss. We add the gradient
$\nabla_{x_{t}}\mathcal{L}_{m}^{\uparrow}$ (or
$\nabla_{x_{t}}\mathcal{L}_{m}^{\downarrow}$, respectively) during the DDIM
update in Eq. 6:

$\bm{x}_{t-1}=\sqrt{\bar{\alpha}_{t-1}}\hat{\bm{x}}_{0|t}^{\prime}+\sqrt{1-\bar{\alpha}_{t-1}-\eta^{2}\tilde{\beta}_{t}^{2}}(\bm{\epsilon}_{\theta}^{(t)}(x_{t})+\gamma\nabla_{\bm{x}_{t}}\mathcal{L}_{m}^{\uparrow/\downarrow})+\eta\tilde{\beta}_{t}\bm{\epsilon}\text{,}$
(13)

where the parameter $\gamma$ makes sure that the $l_{2}$ norm of
$\nabla_{x_{t}}\mathcal{L}_{m}^{\uparrow/\downarrow}$ stays small with respect
to $\bm{\epsilon}_{\theta}$. We set $\gamma$ to the following value:

$\gamma=\begin{cases}b\cdot||\bm{\epsilon}_{\theta}^{(t)}||_{2}/||\nabla_{\bm{x}_{t}}\mathcal{L}_{m}^{\uparrow/\downarrow}||_{2}&\text{if
$||\nabla_{\bm{x}_{t}}\mathcal{L}_{m}^{\uparrow/\downarrow}||_{2}>b||\bm{\epsilon}_{\theta}^{(t)}||_{2}$}\\\
1&\text{otherwise}\text{,}\end{cases}$ (14)

where $b$ is a small constant that we set to $b=0.005$ in our experiments. An
ablation study on the effects of the parameter value can be found in the
Appendix. A summary of our method is shown in Fig. 2.

For the baseline RR, we sample 16 images per slice and select
$x_{0}^{\uparrow}$ and $x_{0}^{\downarrow}$ for every segmentation class as
the reconstruction leading to the highest segmentation volume
($V_{seg}^{\uparrow}$) and lowest segmentation volume
($V_{seg}^{\downarrow}$), respectively.

Figure 4: Analysis of different acceleration factors. For RR, we see that the
uncertainty does not increase, even though very high accelerations are tested.
Our method, in contrast, generates more reliable uncertainty boundaries.

###  3.3 Data and models

We train and test our method on the publicly available SKM-TEA [8] dataset and
focus only on the $E1$-echoes. As the spacing of the SKM-TEA data is very
fine-grained (0.3125 mm in the axial ($=$sampling) direction), neighboring
images tend to be similar. We therefore save computational resources during
testing by only considering every 8th slice in the relevant regions,
corresponding to a spacing of 2.5mm between the slices. In total we analyze
$869$ slices for each acceleration factor and each method.

All experiments are based on the same diffusion model. We train the diffusion
model on the fully-sampled images of the SKM-TEA data. We normalize the data
by dividing by the 0.99 quantile of every slice. The diffusion model is based
on the one introduced in [9]. We train the model for 1 million steps using a
batch-size of 2 on a 3080TI GPU. For segmentation, we train a U-Net on the
fully-sampled MRI-data in axial direction. Due to the strong class-imbalance
on the axial imaging plane, where only a minority of slices show elements of
tibial cartilage or meniscus as these parts tend to be relatively flat, we
perform a foreground oversampling of 40%. To reduce overfitting, we apply
mirroring and Gaussian Noise augmentations during training.

Figure 5: Results. Median score and percentile boundaries (25%, 75%) for
different acceleration factors. Note that the upper and lower bound
segmentation of our method are mostly correctly capturing the inherent
uncertainty even for high acceleration factors: we have high precision (small
amount of false-positives in $S^{\downarrow}$) and recall (small amount of
false-negatives in $S^{\uparrow}$). The ratio of uncertain volume
$V_{unc}=V_{seg}^{\uparrow}-V_{seg}^{\downarrow}$ over $V_{seg}^{\uparrow}$
also increases, better reflecting the true underlying uncertainty.

##  4 Experiments and results

We perform experiments on 4 acceleration factors: $4\times$, $8\times$,
$12\times$, and $16\times$, analyzing the behavior of the uncertainty
boundaries of SGR and RR. We assess the quality of the generated segmentations
$S^{\downarrow}$ and $S^{\uparrow}$ by calculating the precision for
$S^{\downarrow}$, recall for $S^{\uparrow}$ and the uncertainty volume
$V_{unc}$. We focus on precision on $S^{\downarrow}$, as we want to avoid
overestimating the segmentation and creating false-positives, and measure
recall for $S^{\uparrow}$, because we want to avoid segmentations containing
false-negatives. An example image where RR overestimates $S^{\downarrow}$
while our method predicts accurate boundaries is shown in Fig. 3. A similar
observation can be found in Fig. 5, where our method sustains high precision
and recall values, while these scores degrade with higher acceleration for the
baseline method RR. Reconstruction quality as measured by SSIM and PSNR
remains very high for SGR compared to RR, as shown in Tab. 1.

In Fig. 4, we see that the uncertainty boundary increases with higher
acceleration for our method, while it stays approximately constant for the
baseline method RR. This result is also shown in Fig. 5, where the ratio
$V_{unc}/V_{seg}^{\uparrow}$ increases faster for SGR compared to RR. Given
the relatively low precision and recall scores for RR and the higher scores
for our method, one can conclude that RR underestimates the uncertainty
boundaries, while this problem is largely reduced for the uncertainty
boundaries estimated by SGR.

##  5 Conclusion and Discussion

We have introduced a novel method for segmentation-guided MRI reconstruction
and have shown that it can generate meaningfully diverse reconstructions that
lead to upper and lower bound segmentations. These segmentation can offer an
intuitive way of understanding the segmentation uncertainty that is caused by
the ill-posed nature of the reconstruction problem. As our method also offers
more reliable uncertainty boundaries compared to the standard method of
repeated reconstruction, we believe that it can be helpful in clinical
practice, whenever a trustworthy uncertainty estimation is required.

Table 1: Quantitative results for experiments for different accelerations. All reported segmentation metrics are calculated on the 3D-volumes of all 4 segmentation classes. Median values are reported. SSIM and PSNR values are calculated as the average over all reconstructed volumes. Note that it is not clear if larger or smaller ratios of $V_{unc}/V^{\uparrow}$ are preferable, due to the inherent uncertainty of the reconstruction problem. |  | Segmentation | Reconstruction
---|---|---|---
Acc. | Method | Prec. $S^{\downarrow}$ | Recall $S^{\uparrow}$ | $V_{unc}/V^{\uparrow}$ | SSIM | PSNR
4x | SGR (ours) | 0.992 | 0.995 | $0.126$ | $0.914$ | 34.6
RR | 0.982 | 0.989 | 0.045 | 0.915 | 34.6
8x | SGR (ours) | 0.989 | 0.992 | 0.322 | 0.857 | 30.8
RR | 0.971 | 0.976 | 0.075 | 0.857 | 30.9
12x | SGR (ours) | 0.984 | 0.991 | 0.473 | 0.818 | 29.1
RR | 0.961 | 0.966 | 0.096 | 0.820 | 29.1
16x | SGR (ours) | 0.980 | 0.990 | 0.586 | 0.788 | 27.9
RR | 0.951 | 0.957 | 0.113 | 0.789 | 28.0

In future work, we want to analyze if robust segmentation networks introduced
in [7] can help to reduce $V_{unc}$ while maintaining high precision and
recall scores. In this work, however, we have focused on more commonly used
neural networks, to show the easy applicability of our method for existing
pipelines.

####  5.0.1 Acknowledgements

Funded by the Deutsche Forschungsgemeinschaft (DFG, German Research
Foundation) under Germany’s Excellence Strategy – EXC number 2064/1 – Project
number 390727645. The authors thank the International Max Planck Research
School for Intelligent Systems (IMPRS-IS) for supporting Jan Nikolas Morshuis.

####  5.0.2

The authors have no competing interests to declare.

## References

  * [1] Acar, M., Çukur, T., Öksüz, İ.: Segmentation-aware mri reconstruction. In: Machine Learning for Medical Image Reconstruction. pp. 53–61. Springer International Publishing, Cham (2022)
  * [2] Calivá, F., Leynes, A.P., Shah, R., Upadhyay Bharadwaj, U., Majumdar, S., Larson, P.E.Z., Pedoia, V.: Breaking speed limits with simultaneous ultra-fast mri reconstruction and tissue segmentation. In: MIDL. Proceedings of Machine Learning Research, PMLR (06–08 Jul 2020)
  * [3] Cheung, M.Y., Netherton, T.J., Court, L.E., Veeraraghavan, A., Balakrishnan, G.: Metric-guided image reconstruction bounds via conformal prediction. arXiv preprint arXiv:2404.15274 (2024)
  * [4] Chung, H., Lee, S., Ye, J.C.: Decomposed diffusion sampler for accelerating large-scale inverse problems. In: ICLR (2024)
  * [5] Chung, H., Ye, J.C.: Score-based diffusion models for accelerated mri. Medical Image Analysis p. 102479 (2022)
  * [6] Cohen, N., Manor, H., Bahat, Y., Michaeli, T.: From posterior sampling to meaningful diversity in image restoration. In: ICLR (2024)
  * [7] Croce, F., Singh, N.D., Hein, M.: Towards reliable evaluation and fast training of robust semantic segmentation models. In: ECCV (2024)
  * [8] Desai, A.D., Schmidt, A.M., Rubin, E.B., Sandino, C.M., Black, M.S., Mazzoli, V., Stevens, K.J., Boutin, R., Re, C., Gold, G.E., et al.: Skm-tea: A dataset for accelerated mri reconstruction with dense image labels for quantitative clinical evaluation. In: NeurIPS Datasets and Benchmarks Track (Round 2) (2021)
  * [9] Dhariwal, P., Nichol, A.Q.: Diffusion models beat GANs on image synthesis. In: NeurIPS (2021)
  * [10] Efron, B.: Tweedie’s formula and selection bias. Journal of the American Statistical Association 106(496), 1602–1614 (2011)
  * [11] Feiner, L.F., Menten, M.J., Hammernik, K., Hager, P., Huang, W., Rueckert, D., Braren, R.F., Kaissis, G.: Propagation and Attribution of Uncertainty in Medical Imaging Pipelines, pp. 1–11. Springer Nature Switzerland, Cham (2023)
  * [12] Fischer, P., Thomas, K., Baumgartner, C.F.: Uncertainty estimation and propagation in accelerated mri reconstruction. In: Uncertainty for Safe Utilization of Machine Learning in Medical Imaging. pp. 84–94. Springer Nature Switzerland, Cham (2023)
  * [13] Ho, J., Jain, A., Abbeel, P.: Denoising diffusion probabilistic models. NeurIPS (2020)
  * [14] Jalal, A., Arvinte, M., Daras, G., Price, E., Dimakis, A., Tamir, J.: Robust compressed sensing MRI with deep generative priors. In: NeurIPS (2021)
  * [15] Küstner, T., Hammernik, K., Rueckert, D., Hepp, T., Gatidis, S.: Predictive uncertainty in deep learning–based mr image reconstruction using deep ensembles: Evaluation on the fastmri data set. Magnetic Resonance in Medicine (2024)
  * [16] Schlemper, J., Castro, D.C., Bai, W., Qin, C., Oktay, O., Duan, J., Price, A.N., Hajnal, J., Rueckert, D.: Bayesian deep learning for accelerated mr image reconstruction. In: Machine Learning for Medical Image Reconstruction. pp. 64–71. Springer International Publishing, Cham (2018)
  * [17] Song, J., Meng, C., Ermon, S.: Denoising diffusion implicit models. In: ICLR (2021)
  * [18] Song, Y., Shen, L., Xing, L., Ermon, S.: Solving inverse problems in medical imaging with score-based generative models. In: ICLR (2022)
  * [19] Song, Y., Sohl-Dickstein, J., Kingma, D.P., Kumar, A., Ermon, S., Poole, B.: Score-based generative modeling through stochastic differential equations. In: ICLR (2021)
  * [20] Sriram, A., Zbontar, J., Murrell, T., Defazio, A., Zitnick, C.L., Yakubova, N., Knoll, F., Johnson, P.: End-to-end variational networks for accelerated mri reconstruction. In: MICCAI. Springer (2020)
  * [21] Sui, B., Lv, J., Tong, X., Li, Y., Wang, C.: Simultaneous image reconstruction and lesion segmentation in accelerated mri using multitasking learning. Medical Physics 48(11), 7189–7198 (2021)
  * [22] Tolpadi, A.A., Bharadwaj, U., Gao, K.T., Bhattacharjee, R., Gassert, F.G., Luitjens, J., Giesler, P., Morshuis, J.N., Fischer, P., Hein, M., et al.: K2s challenge: From undersampled k-space to automatic segmentation. Bioengineering 10(2), 267 (2023)
  * [23] Wu, Z., Yin, T., Sun, Y., Frost, R., van der Kouwe, A., Dalca, A.V., Bouman, K.L.: Learning task-specific strategies for accelerated mri. IEEE Transactions on Computational Imaging (2024)
  * [24] Xie, Y., Li, Q.: Measurement-conditioned denoising diffusion probabilistic model for under-sampled medical image reconstruction. In: MICCAI. Springer Nature Switzerland, Cham (2022)
  * [25] Zbontar, J., Knoll, F., Sriram, A., Murrell, T., Huang, Z., Muckley, M.J., Defazio, A., Stern, R., Johnson, P., Bruno, M., et al.: fastmri: An open dataset and benchmarks for accelerated mri. arXiv preprint arXiv:1811.08839 (2018)

##  Appendix 0.A Supplementary Material

###  0.A.1 Ablation study for the weighting hyperparameter $b$

Figure 6: Ablation study with 16x acceleration on the hyperparameter $b$,
which is responsible to weight the influence of
$\nabla_{x_{t}}\mathcal{L}_{m}^{\uparrow/\downarrow}$ in the calculation of
$x_{t-1}$ (see Eqs. 13 and 14 in the main paper). Note that this parameter
largely influences the uncertainty volume $V_{unc}$, as well as the Precision
of $S^{\downarrow}$ and the Recall of $S^{\uparrow}$.


# Incoherent diffractive dijet production and gluon Bose enhancement in the
nuclear wave function

Tiyasa Kar<EMAIL_ADDRESS>Department of Physics, North Carolina State
University, Raleigh, NC 27695, USA  Alexander Kovner
<EMAIL_ADDRESS>Physics Department, University of Connecticut, 2152
Hillside Road, Storrs, CT 06269, USA  Ming Li<EMAIL_ADDRESS>Department of
Physics, The Ohio State University, Columbus, OH 43210 , USA  Vladimir V.
Skokov<EMAIL_ADDRESS>Department of Physics, North Carolina State
University, Raleigh, NC 27695, USA

###### Abstract

We investigate the effect of gluon Bose enhancement in the nuclear wave
function on the dijet production in incoherent diffractive processes in DIS
and ultraperipheral collisions. We demonstrate that Bose enhancement leads to
an enhancement of diffractive dijet production cross section when the
transverse momenta of the two jets are aligned at zero relative angle. This
enhancement is maximal when the magnitude of the transverse momenta of the two
jets are equal, and disappears rather quickly as a function of the ratio of
the two momenta. We study both the dilute limit and fully nonlinear dense
regime where the nuclear wave function is evolved with the leading order
JIMWLK equation. In both cases we observe a visible effect, with it being
enhanced by the evolution due to the dynamical generation of the color
neutralization scale.

##  I Introduction

The investigation of dijet production in Deep Inelastic Scattering (DIS) is
one of the key measurements [1] of the upcoming Electron-Ion Collider (EIC),
as it provides insights into the phenomenon of gluon saturation at high
energies [2, 3, 4, 5, 6, 7, 8].

Most of the studies so far, with rare exceptions [9], have been concerned with
the nearly back-to-back dijets [10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20,
21, 22, 23, 24, 25, 26, 9, 27, 28, 29, 30, 31]. Recently the dijet calculation
in DIS in the Color Glass Condensate (CGC) formalism was extended to NLO order
in Refs. [32, 33, 34] with particular attention to understanding the interplay
of high energy evolution and the Sudakov suppression.

In this paper we are interested in a particular aspect of dijet production,
namely can it be used as a probe of Bose enhancement of gluons in the nuclear
wave function. This prompts us to consider both inclusive and diffractive
dijet production [35, 36, 37, 38, 39, 40, 9, 41, 42, 43, 44, 31]; and also to
go beyond the back-to-back kinematics.

In two recent papers [45, 46], we have studied an observable that allows to
directly probe Bose-Einstein correlations between gluons in the nucleus [47,
48, 49]. These correlations are very interesting since they provide nontrivial
information about two particle correlations in the nuclear wave function [50,
51, 52, 53, 54] which may significantly extend our understanding beyond the
conventional single particle distributions traditionally probed in DIS
experiments [18, 37, 24, 9]. The observable we considered was rather
challenging experimentally, as it requires the measurement of three jets and
correlations between them. In this work, we argue that gluon Bose-Einstein
correlations can be measured in a simpler setup in diffractive dijet
production.

Quantum statistics that is responsible for the universal source of
correlations between (identical) gluons is Bose enhancement. In the hadron
wave function this effect leads to the amplification in the number of pairs of
gluons with the same quantum numbers and small relative momentum, see more in
Ref. [47, 48, 49].

The reason we expect dijet production to be sensitive to the Bose correlation
is rather transparent. At high energy, in the infinite momentum frame, the
virtual photon fluctuates into a quark-antiquark pair (dipole) with transverse
momenta ${\underline{q}}_{1}$ and ${\underline{q}}_{2}$, which scatters on the
gluon field of the fast-moving hadron target. In the final state this leads to
two (quark-anti-quark) jets with the transverse momenta ${\underline{p}}_{1}$
and ${\underline{p}}_{2}$. In the hadronic wave function prior to the
scattering gluons are correlated via Bose enhancement. That is, the
probability of finding two gluons is enhanced if they have the same color and
transverse momentum (${\underline{k}}_{1}\approx{\underline{k}}_{2}$). Thus
there is an enhanced probability that as a result of scattering, both the
quark and the antiquark absorb target gluons with the same transverse momentum
${\underline{k}}_{1}$. Since the two gluons are also in the color singlet
state, as is the original dipole one can further isolate this contribution by
considering diffractive dijet production. The two jets may be away from the
back-to-back regime and may even have almost the same transverse momenta
(while longitudinal momentum can be very different), if the exchange with the
target is due to hard gluons.

The dijet observable is much simpler and better experimentally accessible than
the trijet we have considered before. An additional advantage is that in the
diffractive configuration there is no radiation outside of the dijet. Thus the
momentum imbalance of the dijet is not affected by the Sudakov radiation [55,
16] and is determined entirely by the momentum exchange with the target.

In this paper we study numerically the diffractive dijet production cross
section as a function of the relative angle between the transverse momenta of
the two jets in the CGC approach. We start with the dilute regime where the
calculation can be done to lowest order in the density of the target. We then
study the full nonlinear case where the target is evolved by the leading order
JIMWLK evolution.

For perturbative calculation the target averaging is done utilizing the
McLerran-Venugopalan (MV) model [56, 57]. For comparison we also introduce an
IR regulator which enforces color neutralization in the target within a finite
transverse distance $1/m$. The initial condition for JIMWLK evolution is given
by the same MV model. JIMWLK evolution generates the color neutralization
scale dynamically. The inverse neutralization scale $m$ grows with rapidity
proportionally to the saturation momentum.

In both cases we find an enhancement of the production for small relative
angle. This enhancement is most significant when the transverse momenta of the
two jets are equal in magnitude, and disappears when the ratio of the two
momenta is of order 1.5-2. The effect is enhanced when the nonvanishing color
neutralization scale, and is more pronounced for the JIMWLK evolved results.

The paper is structured as follows. In Section II we recap the calculation of
the dijet cross section within the CGC approach. In Section III we review the
Bose enhancement in the CGC wavefunction as it appeared in calculations of
particle production in CGC. In Section IV we calculate the diffractive dijet
production in the dilute limit. In Section V we present the results of the
full nonlinear, JIMWLK evolved cross section. Finally in Section VI we briefly
discuss our results.

##  II Inclusive and incoherent diffractive dijet productions

In this section, we review the differential cross section of DIS dijet
production within the Color Glass Condensate (CGC) formalism [6, 7, 8]. As we
alluded to in the introduction, in this framework, the incident photon
undergoes a fluctuation into a quark-antiquark pair, which then interacts with
the nucleus. We consider the photon to be right-moving while the nuclear
target moves to the left. We set the reference frame such that the photon has
a large longitudinal momentum $p^{+}$ and zero transverse momentum; photon’s
virtuality is denoted by $Q^{2}$. The target nucleus carries a large $P^{-}$
component, that is, we have (neglecting the nucleus mass)

$p^{\mu}=\bigg{(}p^{+},\frac{-Q^{2}}{2p^{+}},{\underline{0}}\bigg{)}$

for the photon and

$P^{\mu}=\bigg{(}0,P^{-},{\underline{0}}\bigg{)}$

for the nucleus. At high energy, it is convenient to work in the so-called
dipole picture of DIS, where the photon splits into quark-anti-quark pair with
the lifetime of the dipole being much longer than the interaction time with
the nucleus. We work in the eikonal approximation, in which the transverse
separation between the quark and antiquark remains unchanged while the dipole
traverses through the nucleus. Corrections to eikonal approximation are
suppressed by powers of collision energy $\sqrt{s}$ [58, 59, 60, 61] and
although this suppression is not overwhelming at EIC we will not consider such
corrections in the present work.

The expressions for the differential dijet production cross-sections can be
found in Ref. [12]:

$\displaystyle\frac{d\sigma^{\gamma_{T,L}^{\ast}A\rightarrow
q\bar{q}X}}{d^{3}k_{1}d^{3}k_{2}}=$ $\displaystyle
N_{c}\alpha_{em}e_{q}^{2}\delta(p^{+}-k_{1}^{+}-k_{2}^{+})\int\frac{\text{d}^{2}{\underline{x}}_{1}}{(2\pi)^{2}}\frac{\text{d}^{2}{\underline{x}}_{1}^{\prime}}{(2\pi)^{2}}\frac{\text{d}^{2}{\underline{x}}_{2}}{(2\pi)^{2}}\frac{\text{d}^{2}{\underline{x}}_{2}^{\prime}}{(2\pi)^{2}}e^{-i{\underline{k}}_{1}\cdot({\underline{x}}_{1}-{\underline{x}}_{1}^{\prime})}e^{-i{\underline{k}}_{2}\cdot({\underline{x}}_{2}-{\underline{x}}_{2}^{\prime})}$
$\displaystyle\times\sum_{\lambda\alpha\beta}\psi^{T,L\lambda}_{\alpha\beta}({\underline{x}}_{1}-{\underline{x}}_{2})\psi^{T,L\lambda*}_{\alpha\beta}({\underline{x}}_{1}^{\prime}-{\underline{x}}_{2}^{\prime})\,\mathcal{N}.$
(1)

The wave functions of the transverse and longitudinal photon are

$\psi^{T\lambda}_{\alpha\beta}(p^{+},z,{\underline{r}})=2\pi\sqrt{\frac{2}{p^{+}}}\begin{cases}i\epsilon_{f}K_{1}(\epsilon_{f}|{\underline{r}}|)\frac{{\underline{r}}\cdot{\underline{\epsilon}}^{(1)}}{|{\underline{r}}|}\left[\delta_{\alpha+}\delta_{\beta+}\bar{z}+\delta_{\alpha-}\delta_{\beta-}z\right],&\text{if
}\lambda=1\\\
i\epsilon_{f}K_{1}(\epsilon_{f}|{\underline{r}}|)\frac{{\underline{r}}\cdot{\underline{\epsilon}}^{(2)}}{|{\underline{r}}|}\left[\delta_{\alpha-}\delta_{\beta-}\bar{z}+\delta_{\alpha+}\delta_{\beta+}z\right],&\text{if
}\lambda=2\end{cases}$ (2)

for the transverse ($T$) and

$\psi^{L\lambda}_{\alpha\beta}(p^{+},z,{\underline{r}})=2\pi\sqrt{\frac{4}{p^{+}}}z\bar{z}QK_{0}(\epsilon_{f}|{\underline{r}}|)\delta_{\alpha\beta}$
(3)

for the longitudinal ($L$) polarizations. The variable $z$ represents the
momentum fraction carried by the quark, and
${\underline{r}}={\underline{x}}_{1}-{\underline{x}}_{2}$ represents the
transverse separation of the quark-antiquark pair. The helicities of the
(massless) quark and antiquark are denoted by $\alpha$ and $\beta$.
Additionally, it is convenient to introduce $\epsilon_{f}^{2}=z\bar{z}Q^{2}$.
The function $K_{n}$ is the modified Bessel function of the second kind of
$n^{\text{th}}$ order. In what follows, we will only consider longitudinal
polarization of the virtual photon.

In Eq.(II), $\mathcal{N}$ is the factor expressing the dipole scattering off
the nucleus. For inclusive dijet production it is

$\mathcal{N}_{\text{inclusive}}=1+S^{(4)}_{x_{g}}({\underline{x}}_{1},{\underline{x}}_{2};{\underline{x}}_{2}^{\prime},{\underline{x}}_{1}^{\prime})-S^{(2)}_{x_{g}}({\underline{x}}_{1},{\underline{x}}_{2})-S^{(2)}_{x_{g}}({\underline{x}}_{2}^{\prime},{\underline{x}}_{1}^{\prime}),$
(4)

while for incoherent diffractive dijet production,

$\mathcal{N}_{\text{incoherent
diffractive}}=S^{(2,2)}_{x_{g}}({\underline{x}}_{1},{\underline{x}}_{2};{\underline{x}}_{1}^{\prime},{\underline{x}}_{2}^{\prime})-S^{(2)}_{x_{g}}({\underline{x}}_{1},{\underline{x}}_{2})S^{(2)}_{x_{g}}({\underline{x}}_{1}^{\prime},{\underline{x}}_{2}^{\prime}).$
(5)

The quadrupole and dipole terms arising in ${\cal N}$ of Eq. (II),
respectively, are as follows:

$\displaystyle
S^{(4)}_{x_{g}}({\underline{x}}_{1},{\underline{x}}_{2};{\underline{x}}_{2}^{\prime},{\underline{x}}_{1}^{\prime})$
$\displaystyle=\frac{1}{N_{c}}{\rm Tr}\langle
V^{\dagger}({\underline{x}}_{2})V({\underline{x}}_{1})\left[V^{\dagger}({\underline{x}}^{\prime}_{2})V({\underline{x}}^{\prime}_{1})\right]^{\dagger}\rangle\,,$
(6) $\displaystyle S^{(2)}_{x_{g}}({\underline{x}}_{1},{\underline{x}}_{2})$
$\displaystyle=\frac{1}{N_{c}}{\rm Tr}\langle
V^{\dagger}({\underline{x}}_{2})V({\underline{x}}_{1})\rangle$ (7)

where the Wilson lines $V({\underline{x}})$ (and
$V^{\dagger}({\underline{x}})$) encode multiple scattering of the quark (anti-
quark) off the gluon field of the nucleus. Here $\langle\ldots\rangle$ denotes
averaging over the valence target ensemble.

The additional quantity that appears in the incoherent diffractive cross
section is

$S^{(2,2)}_{x_{g}}({\underline{x}}_{1},{\underline{x}}_{2};{\underline{x}}_{1}^{\prime},{\underline{x}}_{2}^{\prime})=\frac{1}{N_{c}^{2}}\left\langle{\rm
Tr}\left[V^{\dagger}({\underline{x}}_{2})V({\underline{x}}_{1})\right]{\rm
Tr}\left[V^{\dagger}({\underline{x}}^{\prime}_{1})V({\underline{x}}^{\prime}_{2})\right]\right\rangle\,.$
(8)

##  III Bose enhancement

###  III.1 Soft gluon density matrix and averaging over target configurations

CGC is an effective field theory based on the separation of the degrees of
freedom into valence and soft sectors as quantified by the fraction of the
longitudinal momentum of the hadron they carry. The valence charges are
treated as static sources of faster soft degrees of freedom.

The conventional approach in CGC is to first average over the soft sector and
then over the valence target configurations – equations in the previous
section were obtained in this precise way. Although this approach is most
convenient for practical calculations it is not very illuminating as far as
the quantum correlations in the target are concerned. To analyze the latter,
one, first, averages over valence sector and only then considers expectations
values of operators of interest in the soft sector.

The CGC wave function [62, 63] is given by

$|\psi\rangle=|s_{v}\rangle\otimes|v\rangle\,,$ (9)

where $|v\rangle$ is the state vector characterizing the valence degrees of
freedom and $|s_{v}\rangle$ is the vacuum of the soft fields in the presence
of the valence source.

In the leading perturbative order

$|s_{v}\rangle=\mathcal{C}|0\rangle,\quad{\cal C}=\exp\left\\{2i{\rm
tr}\int_{{\underline{k}}}b^{i}({\underline{k}})\phi^{a}_{i}({\underline{k}})\right\\}\,$
(10)

where

$\phi_{i}({\underline{k}})\equiv
a_{i}^{+}({\underline{k}})+a_{i}(-{\underline{k}})\,$ (11)

is the soft gluon field operator with $a^{+}_{i}$, $a_{i}$ being the gluon
creation and annihilation operators. The background field $b^{i}_{a}$ is
determined by the valence color charge density $\rho$ via:

$b^{i}_{a}({\underline{k}})=g\rho_{a}({\underline{k}})\frac{i{\underline{k}}_{i}}{k^{2}}+{\cal
O}(\rho^{2})\,.$ (12)

In this and the next section, for the sake of simplicity of presentation, we
consider the dilute approximation in the target, that is the leading
contribution of the color charge density. At the leading order in
$\rho({\underline{k}})$, only gluons with the longitudinal polarization
contribute to ${\cal C}$ and $|s_{v}\rangle$.

The valence wave function $|v\rangle$ is customarily modeled in the so-called
McLerran-Venugopalan (MV) model as [57, 56]. The MV model can be formulated in
terms of diagonal matrix elements of the valence density matrix in the color
charge density basis:

$\langle\rho|v\rangle\langle
v|\rho\rangle=Ne^{-\int_{{\underline{k}}}\frac{1}{2\mu^{2}}\rho_{a}({\underline{k}})\rho^{*}_{a}({\underline{k}})}\,,$
(13)

where $N$ is the normalization factor and the parameter $\mu^{2}$ determines
the average color charge density in the valence wave function.

The full hadron density matrix reads

$\hat{\rho}=|v\rangle\otimes|s_{v}\rangle\langle s_{v}|\otimes\langle v|\,.$
(14)

Integrating out the valence modes yields the reduced density matrix for the
soft sector

$\hat{\rho}_{r}={\rm Tr}_{\rho}\hat{\rho}\equiv\int
D\rho\,\langle\rho|\hat{\rho}|\rho\rangle=\int
D\rho\,\langle\rho|v\rangle\,|s_{v}\rangle\langle s_{v}|\,\langle
v|\rho\rangle\,.$ (15)

In the particle number representation it was derived in Ref. [63]. The
nonvanishing matrix elements of $\hat{\rho}_{r}$ between states in the
momentum space Fock basis are

$\displaystyle\rho_{n,m,\alpha,\beta}\equiv\langle
n_{c}({\underline{q}}),m_{c}(-{\underline{q}})|\hat{\rho}_{r}({\underline{q}})|\alpha_{c}({\underline{q}}),\beta_{c}(-{\underline{q}})\rangle=(1-R)\frac{\left(n+\beta\right)!}{\sqrt{n!m!\alpha!\beta!}}\left(\frac{R}{2}\right)^{n+\beta}\delta_{\left(n+\beta\right),\left(m+\alpha\right)}\,,$
(16) $\displaystyle
R=\left(1+\frac{{\underline{q}}^{2}}{2g^{2}\mu^{2}}\right)^{-1}$ (17)

with other elements being zero. The matrix elements are explicitly diagonal in
color and momentum. The state
$|n_{c}({\underline{q}}),m_{c}(-{\underline{q}})\rangle$ denotes $n$ gluons
with momentum ${\underline{q}}$ and $m$ gluons with momentum
$-{\underline{q}}$.

###  III.2 Bose enhancement in the soft density matrix

Consider the correlator (only longitudinal polarization of gluons is included)

$\displaystyle D({\underline{k}},{\underline{p}})={\rm
Tr}\left(\hat{\rho}_{r}a_{b}^{+}({\underline{k}})a_{c}^{+}({\underline{p}})a_{b}({\underline{k}})a_{c}({\underline{p}})\right).$
(18)

This measures the correlation between two gluons drawn from the target.

This can be easily calculated in the MV model. Consider first

$\displaystyle\langle
a_{a}^{+}({\underline{k}}_{1})a_{b}({\underline{k}}_{2})\rangle={\rm
Tr}\left(\hat{\rho}_{r}a_{a}^{+}({\underline{k}}_{1})a_{b}({\underline{k}}_{2})\right)=(2\pi)^{2}\delta^{(2)}({\underline{k}}_{1}-{\underline{k}}_{2})\delta_{ab}\sum_{n,m}n\rho_{n,m,n,m}=(2\pi)^{2}\delta^{(2)}({\underline{k}}_{1}-{\underline{k}}_{2})\delta_{ab}\frac{g^{2}\bar{\mu}^{2}}{k_{1}^{2}}\,.$
(19)

The sum is computed in Appendix A.

Additionally

$\displaystyle\langle
a_{a}^{+}({\underline{k}}_{1})a^{+}_{b}({\underline{k}}_{2})\rangle=(2\pi)^{2}\delta^{(2)}({\underline{k}}_{1}-{\underline{k}}_{2})\delta_{ab}\sum_{n,m}\sqrt{n+1}\sqrt{m+1}\rho_{n,m,n+1,m+1}=(2\pi)^{2}\delta^{(2)}({\underline{k}}_{1}+{\underline{k}}_{2})\delta_{ab}\frac{g^{2}\bar{\mu}^{2}}{k_{1}^{2}}\,$
(20)

with the same result for $\langle
a_{a}({\underline{k}}_{1})a_{b}({\underline{k}}_{2})\rangle$.

In Ref. [64], it was shown that the density matrix $\rho_{r}$ is Gaussian,
that is $\hat{\rho}_{r}=\exp\left[-\int
d^{2}q\,\beta\omega({\underline{q}})c^{+}({\underline{q}})c({\underline{q}})\right]$
where $c({\underline{q}})$ is a linear combination of $a^{+}({\underline{q}})$
and $a(-{\underline{q}})$. Therefore we can use Wick theorem to find the
expectation value of $D({\underline{k}},{\underline{p}})$. We should keep in
mind that in addition to the standard Wick contraction $\langle
a^{\dagger}a\rangle$, also the contractions $\langle
a^{\dagger}a^{\dagger}\rangle$ and $\langle aa\rangle$ are nonvanishing and
are given by (20). We get

$\displaystyle D({\underline{k}},{\underline{p}})$ $\displaystyle=\left\langle
a_{b}^{+}({\underline{k}})a_{c}^{+}({\underline{p}})a_{b}({\underline{k}})a_{c}({\underline{p}})\right\rangle$
$\displaystyle=\left\langle
a_{b}^{+}({\underline{k}})a_{b}({\underline{k}})\right\rangle\left\langle
a_{c}^{+}({\underline{p}})a_{c}({\underline{p}})\right\rangle+\left\langle
a_{b}^{+}({\underline{k}})a_{c}^{+}({\underline{p}})\right\rangle\left\langle
a_{b}({\underline{k}})a_{c}({\underline{p}})\right\rangle+\left\langle
a_{b}^{+}({\underline{k}})a_{c}({\underline{p}})\right\rangle\left\langle
a_{c}^{+}({\underline{p}})a_{b}({\underline{k}})\right\rangle\,.$ (21)

Here, as we are interested in density-enhanced terms, we can ignore the order
of the operators, that is we can treat the creation annihilation operators as
commuting. We thus obtain

$\displaystyle D({\underline{k}},{\underline{p}})$
$\displaystyle=\underbrace{\left(S_{\perp}(N_{c}^{2}-1)\frac{g^{2}\bar{\mu}^{2}}{{\underline{k}}^{2}}\right)}_{n({\underline{k}})}\underbrace{\left(S_{\perp}(N_{c}^{2}-1)\frac{g^{2}\bar{\mu}^{2}}{{\underline{p}}^{2}}\right)}_{n({\underline{p}})}$
$\displaystyle+(2\pi)^{2}(N_{c}^{2}-1)S_{\perp}\left(\frac{g^{2}\bar{\mu}^{2}}{{\underline{k}}^{2}}\right)^{2}\left[\delta^{(2)}({\underline{k}}+{\underline{p}})+\delta^{(2)}({\underline{k}}-{\underline{p}})\right]\,.$
(22)

The first term is the “classical” independent and factorizing contribution,
while the last term represents the non-trivial effect of Bose correlations.
Note that the correlation contains the sum of two delta functions. The first
one leads to the enhancement of gluon correlations for back-to-back
kinematics. Although this term is important, in physical processes, it is
often overpowered by the stronger effect of the transverse momentum
conservation. The later term leads to the enhancement of gluons at the same
momentum and this is the effect we are interested in. Note that the
sensitivity to the Bose enhancement starts from the observables quartic in
creation/annihilation operators or, in other words, from the order of
$\bar{\mu}^{4}$.

##  IV Dilute approximation

We now discuss averaging over the target in dilute approximation. The basic
objects we have to average are products of Wilson lines. The Wilson line in
the fundamental representation is

$\displaystyle V({\underline{x}})={\cal
P}\exp\left(ig\int_{-\infty}^{+\infty}dx^{-}t^{a}A_{a}^{+}(x^{-},x^{+}=0,{\underline{x}})\right)\,.$
(23)

In the covariant gauge

$\displaystyle\partial^{2}A_{a}^{+}(x^{-},{\underline{x}})=g\rho_{a}(x^{-},{\underline{x}})\,.$
(24)

which yields

$\displaystyle
A_{a}^{+}(x^{-},{\underline{x}})=g\frac{1}{\partial^{2}}\rho(x^{-},{\underline{x}})=-\frac{g}{2\pi}\int
d^{2}{\underline{y}}\ \ln(|{\underline{x}}-{\underline{y}}|\Lambda)\
\rho(x^{-},{\underline{y}})=\frac{g}{2\pi}\int\frac{d^{2}{\underline{k}}}{(2\pi)^{2}}\frac{1}{k^{2}}\rho(x^{-},{\underline{k}})\,.$
(25)

In the MV model, the correlators of $\rho$ are determined by a Gaussian
distribution with

$\displaystyle\langle\rho^{a}(x^{-},{\underline{x}})\rho^{b}(y^{-},{\underline{y}})\rangle=\delta^{ab}\mu^{2}(x^{-})\delta(x^{-}-y^{-})\delta^{(2)}({\underline{x}}-{\underline{y}})\,.$
(26)

For the field correlators this translates into

$\displaystyle\langle
A^{a}(x^{-},{\underline{x}})A^{b}(y^{-},{\underline{y}})\rangle=\delta^{ab}g^{2}\mu^{2}(x^{-})\delta(x^{-}-y^{-})L({\underline{x}}-{\underline{y}})$
(27)

where

$L({\underline{x}}-{\underline{y}})=\frac{g^{2}}{(2\pi)^{2}}\int
d^{2}z\ln(|{\underline{x}}-{\underline{z}}|\Lambda)\ln(|{\underline{z}}-{\underline{y}}|\Lambda)\,.$

In the dilute limit one expands in the number of gluon exchanges with the
target, which necessitates expansion of the Wilson line. As is clear from the
previous discussion, the lowest order in which the effect of Bose enhancement
is visible is the fourth order in $\mu$, which means that in the expansion of
a single Wilson line we need to keep terms up to fourth order in the field

$\displaystyle V(x)\approx 1$
$\displaystyle+ig\int_{-\infty}^{+\infty}dx^{-}t^{a}A_{a}^{+}(x^{-},{\underline{x}})+(ig)^{2}\int_{-\infty}^{+\infty}dx_{0}^{-}\int_{-\infty}^{x_{0}^{-}}dx_{1}^{-}t^{a}t^{b}A_{a}^{+}(x_{0}^{-},{\underline{x}})A_{b}^{+}(x_{1}^{-},{\underline{x}})$
$\displaystyle+(ig)^{3}\int_{-\infty}^{+\infty}dx_{0}^{-}\int_{-\infty}^{x_{0}^{-}}dx_{1}^{-}\int_{-\infty}^{x_{1}^{-}}dx_{2}^{-}t^{a}t^{b}t^{c}A_{a}^{+}(x_{0}^{-},{\underline{x}})A_{b}^{+}(x_{1}^{-},{\underline{x}})A_{c}^{+}(x_{2}^{-},{\underline{x}})$
$\displaystyle+(ig)^{4}\int_{-\infty}^{+\infty}dx_{0}^{-}\int_{-\infty}^{x_{0}^{-}}dx_{1}^{-}\int_{-\infty}^{x_{1}^{-}}dx_{2}^{-}\int_{-\infty}^{x_{2}^{-}}dx_{3}^{-}t^{a}t^{b}t^{c}t^{d}A_{a}^{+}(x_{0}^{-},{\underline{x}})A_{b}^{+}(x_{1}^{-},{\underline{x}})A_{c}^{+}(x_{2}^{-},{\underline{x}})A_{d}^{+}(x_{3}^{-},{\underline{x}})\,.$
(28)

To make our notation compact, we introduce the averaged quantities,

$\displaystyle\bar{\mu}^{2}=\int_{-\infty}^{+\infty}dx^{-}\mu^{2}(x^{-})$ (29)

and

$\displaystyle\alpha_{a}({\underline{x}})=\int_{-\infty}^{+\infty}dx_{0}^{-}A_{a}^{+}(x_{0}^{-},{\underline{x}})\,.$
(30)

We also notice that the cubic order in the field can only appear in the
expectation values when multiplied by the linear order. We thus can contract
two remaining fields already in the expansion of the Wilson line, that is

$\displaystyle\int_{-\infty}^{+\infty}dx_{0}^{-}\int_{-\infty}^{x_{0}^{-}}dx_{1}^{-}\int_{-\infty}^{x_{1}^{-}}dx_{2}^{-}t^{a}t^{b}t^{c}A_{a}^{+}(x_{0}^{-},{\underline{x}})A_{b}^{+}(x_{1}^{-},{\underline{x}})A_{c}^{+}(x_{2}^{-},{\underline{x}})\to\int_{-\infty}^{+\infty}dx_{0}^{-}\int_{-\infty}^{x_{0}^{-}}dx_{1}^{-}\int_{-\infty}^{x_{1}^{-}}dx_{2}^{-}t^{a}t^{b}t^{c}$
$\displaystyle\quad\times\left(A_{a}^{+}(x_{0}^{-},{\underline{x}})\langle
A_{b}^{+}(x_{1}^{-},{\underline{x}})A_{c}^{+}(x_{2}^{-},{\underline{x}})\rangle+\langle
A_{a}^{+}(x_{0}^{-},{\underline{x}})A_{b}^{+}(x_{1}^{-},{\underline{x}})\rangle
A_{c}^{+}(x_{2}^{-},{\underline{x}})\right.$
$\displaystyle\left.\quad\quad\quad+A_{b}^{+}(x_{1}^{-},{\underline{x}})\langle
A_{a}^{+}(x_{0}^{-},{\underline{x}})A_{c}^{+}(x_{2}^{-},{\underline{x}})\rangle\right)$
$\displaystyle=\frac{C_{f}L({\underline{0}})}{2}\int_{-\infty}^{+\infty}dx_{0}^{-}\int_{-\infty}^{x_{0}^{-}}dx_{1}^{-}t^{a}A_{a}^{+}(x_{0}^{-},{\underline{x}})\mu^{2}(x_{1}^{-})+\frac{C_{f}L({\underline{0}})}{2}\int_{-\infty}^{+\infty}dx_{0}^{-}\int_{-\infty}^{x_{0}^{-}}dx_{2}^{-}t^{a}A_{a}^{+}(x_{2}^{-},{\underline{x}})\mu^{2}(x_{0}^{-})$
$\displaystyle=\frac{C_{f}L({\underline{0}})}{2}\int_{-\infty}^{+\infty}dx_{0}^{-}\int_{-\infty}^{+\infty}dx_{1}^{-}t^{a}A_{a}^{+}(x_{0}^{-},{\underline{x}})\mu^{2}(x_{1}^{-})=\frac{C_{f}\bar{\mu}^{2}L({\underline{0}})}{2}t^{a}\alpha_{a}({\underline{x}})\,.$
(31)

Similar logic allows us to simplify the quartic term (here only one
contraction survives due to color structure):

$\displaystyle\int_{-\infty}^{+\infty}dx_{0}^{-}\int_{-\infty}^{x_{0}^{-}}dx_{1}^{-}\int_{-\infty}^{x_{1}^{-}}dx_{2}^{-}\int_{-\infty}^{x_{2}^{-}}dx_{3}^{-}t^{a}t^{b}t^{c}t^{d}A_{a}^{+}(x_{0}^{-},{\underline{x}})A_{b}^{+}(x_{1}^{-},{\underline{x}})A_{c}^{+}(x_{2}^{-},{\underline{x}})A_{d}^{+}(x_{3}^{-},{\underline{x}})$
$\displaystyle\to$
$\displaystyle\int_{-\infty}^{+\infty}dx_{0}^{-}\int_{-\infty}^{x_{0}^{-}}dx_{1}^{-}\int_{-\infty}^{x_{1}^{-}}dx_{2}^{-}\int_{-\infty}^{x_{2}^{-}}dx_{3}^{-}t^{a}t^{b}t^{c}t^{d}\langle
A_{a}^{+}(x_{0}^{-},{\underline{x}})A_{b}^{+}(x_{1}^{-},{\underline{x}})\rangle\langle
A_{c}^{+}(x_{2}^{-},{\underline{x}})A_{d}^{+}(x_{3}^{-},{\underline{x}})\rangle$
$\displaystyle=$
$\displaystyle\frac{1}{2}\left(\frac{C_{f}\bar{\mu}^{2}L({\underline{0}})}{2}\right)^{2}\,.$
(32)

We thus effectively have

$\displaystyle V({\underline{x}})\approx$
$\displaystyle\,1+\frac{(ig)^{4}}{2}\left(\frac{C_{f}g^{2}\bar{\mu}^{2}L({\underline{0}})}{2}\right)^{2}+igt^{a}\alpha_{a}({\underline{x}})\left(1+(ig)^{2}\frac{C_{f}\bar{g}^{2}\mu^{2}L({\underline{0}})}{2}\right)$
$\displaystyle+(ig)^{2}\int_{-\infty}^{+\infty}dx_{0}^{-}\int_{-\infty}^{x_{0}^{-}}dx_{1}^{-}t^{a}t^{b}A_{a}^{+}(x_{0}^{-},{\underline{x}})A_{b}^{+}(x_{1}^{-},{\underline{x}})\,.$
(33)

Now, to evaluate the diffractive cross-section, it is also useful to compute
the dipole operator to $\mu^{4}$ order

$\displaystyle\frac{1}{N_{c}}{\rm
Tr}\,V^{\dagger}({\underline{y}})V({\underline{x}})$
$\displaystyle=1+(ig)^{4}\left(\frac{C_{f}g^{2}\bar{\mu}^{2}L({\underline{0}})}{2}\right)^{2}$
$\displaystyle-(ig)^{2}\frac{1}{N_{c}}{\rm
Tr}(t_{a}t_{b})\alpha_{a}({\underline{x}})\alpha_{b}({\underline{y}})$
$\displaystyle+(ig)^{2}\frac{1}{N_{c}}{\rm
Tr}(t_{a}t_{b})\int_{-\infty}^{+\infty}dy_{0}^{-}\int^{-\infty}_{y_{0}^{-}}dy_{1}^{-}A_{a}^{+}(y_{0}^{-},{\underline{y}})A_{b}^{+}(y_{1}^{-},{\underline{y}})$
$\displaystyle+(ig)^{2}\frac{1}{N_{c}}{\rm
Tr}(t_{a}t_{b})\int_{-\infty}^{+\infty}dx_{0}^{-}\int_{-\infty}^{x_{0}^{-}}dx_{1}^{-}A_{a}^{+}(x_{0}^{-},{\underline{x}})A_{b}^{+}(x_{1}^{-},{\underline{x}})$
$\displaystyle-(ig)^{4}\frac{1}{N_{c}}{\rm
Tr}(t_{a}t_{b})\left(C_{f}g^{2}\bar{\mu}^{2}L({\underline{0}})\right)\alpha_{a}({\underline{x}})\alpha_{b}({\underline{y}})$
$\displaystyle+(ig)^{3}\frac{1}{N_{c}}{\rm
Tr}\,(t^{b}t^{c}t^{a})\alpha_{a}({\underline{x}})\int_{-\infty}^{+\infty}dy_{0}^{-}\int^{-\infty}_{y_{0}^{-}}dy_{1}^{-}A_{b}^{+}(y_{0}^{-},{\underline{y}})A_{c}^{+}(y_{1}^{-},{\underline{y}})$
$\displaystyle-(ig)^{3}\frac{1}{N_{c}}{\rm
Tr}\,(t^{a}t^{b}t^{c})\alpha_{a}({\underline{y}})\int_{-\infty}^{+\infty}dx_{0}^{-}\int_{-\infty}^{x_{0}^{-}}dx_{1}^{-}A_{b}^{+}(x_{0}^{-},{\underline{x}})A_{c}^{+}(x_{1}^{-},{\underline{x}})$
$\displaystyle+(ig)^{4}\frac{1}{N_{c}}{\rm
Tr}\,(t^{a}t^{b}t^{c}t^{d})\int_{-\infty}^{+\infty}dy_{0}^{-}\int^{-\infty}_{y_{0}^{-}}dy_{1}^{-}A_{a}^{+}(y_{0}^{-},{\underline{y}})A_{b}^{+}(y_{1}^{-},{\underline{y}})$
$\displaystyle\quad\quad\quad\quad\quad\quad\times\int_{-\infty}^{+\infty}dx_{0}^{-}\int_{-\infty}^{x_{0}^{-}}dx_{1}^{-}A_{c}^{+}(x_{0}^{-},{\underline{x}})A_{d}^{+}(x_{1}^{-},{\underline{x}})\,.$
(34)

The last term can be immediately simplified as we did for a single Wilson
line. The cubic term in $A$ involves a trace of three Gell-Mann matrices and,
to this order, will necessarily involve one internal contraction leading to
zero contribution. We thus can safely drop the cubic terms. After this
simplification, we obtain

$\displaystyle\frac{1}{N_{c}}{\rm
Tr}\,V^{\dagger}({\underline{y}})V({\underline{x}})$
$\displaystyle=1-(ig)^{2}\frac{1}{N_{c}}{\rm
Tr}(t_{a}t_{b})\alpha_{a}({\underline{x}})\alpha_{b}({\underline{y}})$
$\displaystyle+(ig)^{2}\frac{1}{N_{c}}{\rm
Tr}(t_{a}t_{b})\int_{-\infty}^{+\infty}dy_{0}^{-}\int^{-\infty}_{y_{0}^{-}}dy_{1}^{-}A_{a}^{+}(y_{0}^{-},{\underline{y}})A_{b}^{+}(y_{1}^{-},{\underline{y}})$
$\displaystyle+(ig)^{2}\frac{1}{N_{c}}{\rm
Tr}(t_{a}t_{b})\int_{-\infty}^{+\infty}dx_{0}^{-}\int_{-\infty}^{x_{0}^{-}}dx_{1}^{-}A_{a}^{+}(x_{0}^{-},{\underline{x}})A_{b}^{+}(x_{1}^{-},{\underline{x}})$
$\displaystyle+\frac{1}{2}\,(ig)^{4}\left(C_{f}g^{2}\bar{\mu}^{2}L({\underline{0}})\right)^{2}-(ig)^{4}\left(C_{f}g^{2}\bar{\mu}^{2}\right)^{2}L({\underline{0}})L({\underline{x}}-{\underline{y}})$
$\displaystyle+\frac{1}{2}(ig)^{4}(C_{f}g^{2}\bar{\mu}^{2})^{2}L^{2}({\underline{x}}-{\underline{y}})$
$\displaystyle=1+\frac{(ig)^{2}}{2}\frac{1}{N_{c}}{\rm
Tr}(t_{a}t_{b})(\alpha_{a}({\underline{y}})-\alpha_{a}({\underline{x}}))(\alpha_{b}({\underline{y}})-\alpha_{b}({\underline{x}}))$
$\displaystyle+\frac{(ig)^{4}(C_{f}g^{2}\bar{\mu}^{2})^{2}}{2}[L({\underline{0}})-L({\underline{x}}-{\underline{y}})]^{2}\,.$
(35)

The diffractive cross-section envolves only a bilinear of the fluctuation of
the dipole amplitude, $\frac{1}{N_{c}}{\rm
Tr}\,V^{\dagger}({\underline{y}})V({\underline{x}})-\left\langle\frac{1}{N_{c}}{\rm
Tr}\,V^{\dagger}({\underline{y}})V({\underline{x}})\right\rangle$. Thus when
computing diffractive production we can drop the first and the last term in
Eq. (IV). For the leading contribution we get simply

$\displaystyle\mathcal{N}_{\text{diffractive}}\approx$
$\displaystyle\frac{g^{4}}{16}\langle(\alpha_{a}({\underline{x}}_{2})-\alpha_{a}({\underline{x}}_{1}))(\alpha_{a}({\underline{x}}_{2})-\alpha_{a}({\underline{x}}_{1}))(\alpha_{b}({\underline{x}}_{2}^{\prime})-\alpha_{b}({\underline{x}}_{1}^{\prime}))(\alpha_{b}({\underline{x}}_{2}^{\prime})-\alpha_{b}({\underline{x}}_{1}^{\prime}))\rangle$
$\displaystyle-\frac{g^{4}}{16}\langle(\alpha_{a}({\underline{x}}_{2})-\alpha_{a}({\underline{x}}_{1}))(\alpha_{a}({\underline{x}}_{2})-\alpha_{a}({\underline{x}}_{1}))\rangle\langle(\alpha_{b}({\underline{x}}_{2}^{\prime})-\alpha_{b}({\underline{x}}_{1}^{\prime}))(\alpha_{b}({\underline{x}}_{2}^{\prime})-\alpha_{b}({\underline{x}}_{1}^{\prime}))\rangle$
(36)

or, evaluating, the averages (see fig.1 for diagrammatic representation)

$\displaystyle\mathcal{N}_{\text{diffractive}}\approx\frac{C_{f}g^{8}\bar{\mu}^{4}}{4N_{c}}(L({\underline{x}}_{1}-{\underline{x}}_{1}^{\prime})-L({\underline{x}}_{1}-{\underline{x}}_{2}^{\prime})-L({\underline{x}}_{1}^{\prime}-{\underline{x}}_{2})+L({\underline{x}}_{2}-{\underline{x}}_{2}^{\prime}))^{2}\,.$
(37)

Similar Wilson line structures in the dilute regime have been analyzed in Ref.
[65].

The cross-section therefore is

$\displaystyle E_{1}E_{2}\frac{d\sigma^{\gamma_{L}^{\ast}A\rightarrow
q\bar{q}X}}{d^{3}k_{1}d^{3}k_{2}}\bigg{|}_{D}=$
$\displaystyle\alpha_{em}e_{q}^{2}Q^{2}z^{2}\bar{z}^{2}\frac{C_{f}g^{8}\bar{\mu}^{4}S_{\perp}}{(2\pi)^{4}}\int\frac{\text{d}^{2}{\underline{q}}}{(2\pi)^{2}}L({\underline{q}})L({\underline{k}}_{1}+{\underline{k}}_{2}-{\underline{q}})$
$\displaystyle\times\Bigg{(}\frac{1}{\epsilon_{f}^{2}+({\underline{k}}_{1}-{\underline{q}})^{2}}-\frac{1}{\epsilon_{f}^{2}+{\underline{k}}_{1}^{2}}+\frac{1}{\epsilon_{f}^{2}+({\underline{k}}_{2}-{\underline{q}})^{2}}-\frac{1}{\epsilon_{f}^{2}+{\underline{k}}_{2}^{2}}\Bigg{)}^{2}\,.$
(38)

$x_{1}\quad x_{1}^{\prime}$$x_{2}\quad x_{2}^{\prime}$ (a) Independent
scattering off the target. This contribution is proportional to
$L({\underline{x}}_{1}-{\underline{x}}^{\prime}_{1})L({\underline{x}}_{2}-{\underline{x}}^{\prime}_{2})$
of Eq. (37).

$x_{1}\quad x_{1}^{\prime}$$\phantom{x_{2}}\quad x_{2}^{\prime}$ (b) Single-
gluon-correlated scattering proportional to
$-L({\underline{x}}_{1}-{\underline{x}}^{\prime}_{2})L({\underline{x}}_{1}-{\underline{x}}^{\prime}_{1})$;
this term somewhat suppresses the back-to-back peak.

$x_{1}\quad x_{1}^{\prime}$$x_{2}\quad x_{2}^{\prime}$ (c) Two gluon
correlated scattering
$L({\underline{x}}_{1}-{\underline{x}}^{\prime}_{2})L({\underline{x}}_{2}-{\underline{x}}^{\prime}_{1})$
with gluons having equal transverse momenta. This contribution causes Bose
Enchancement.

$x_{1}\quad x_{1}^{\prime}$$x_{2}\quad\phantom{x_{2}^{\prime}}$ (d) Two gluon
correlated scattering
$L^{2}({\underline{x}}_{2}-{\underline{x}}^{\prime}_{1})$ with gluons having
opposite transverse momenta. This contribution causes additional enchancemnt
of the back-to-back peak.

Figure 1: Sample of representative diagrams contributing to diffractive dijet
production; note that the final state dijet is required to be a color singlet.

In order to analyze the angular correlations we consider the correlation
function defined as

$\displaystyle
C(\Delta\phi,|{\underline{k}}_{1}|,|{\underline{k}}_{2}|)=\frac{1}{C_{\rm
norm}}\int
d\phi_{1}d\phi_{2}\delta(\Delta\phi-\phi_{1}+\phi_{2})\frac{d\sigma^{\gamma_{L}^{\ast}A\rightarrow
q\bar{q}X}}{d^{3}k_{1}d^{3}k_{2}}\bigg{|}_{D}$ (39)

at a constant value of $z$. The momentum-dependent normalization constant
$C_{\rm norm}$ is fixed by requiring

$\int_{-\pi/2}^{\pi/2}d\Delta\phi\,C(\Delta\phi,|{\underline{k}}_{1}|,|{\underline{k}}_{2}|)=1.$

(a) MV model with no color neutralization, $m\to 0$.

(b) MV model with the color neutralization scale $m=\tilde{Q}_{s}$.

Figure 2: The correlation function (II) for the incoherent diffractive dijet
for $k_{1}=7\tilde{Q}_{s}$ and $\epsilon_{f}=2\tilde{Q}_{s}$.

The results are plotted in Fig. 2 as a function of the relative angle between
the vectors ${\underline{k}}_{1}$ and ${\underline{k}}_{2}$ and their
magnitude. The target averaging was performed with two different models: the
original MV model, and an MV-like model modified to include a finite color
neutralization scale $m$ (see discussion below). The latter is defined by the
field-field correlator modified in the IR:

$\displaystyle
L({\underline{q}})=\frac{1}{{\underline{q}}^{2}({\underline{q}}^{2}+m^{2})}\,.$
(40)

The “mass” $m$ plays the role of the color neutralization scale modeled
through the color charge correlator (see Refs. [66, 67, 45])

$\displaystyle\langle\rho^{a}(x^{-},{\underline{k}})\rho^{b}(y^{-},{\underline{k}}^{\prime})\rangle\propto\frac{k^{2}}{k^{2}+m^{2}}\delta^{ab}\delta(x^{-}-y^{-})\delta({\underline{k}}+{\underline{k}}^{\prime})$
(41)

which ensures that the correlator vanishes at momentum $k\ll m$ (or transverse
separations greater than $1/m$) as ${\underline{k}}^{2}$ (color
neutralization) and that it approaches constant at larger $|{\underline{k}}|$
(a more realistic model should include anomalous dimension at larger $k$).
Equation (40) follows from Eq. (41). Although the original MV model does not
possess the color neutralization property, when evolved with the JIMWLK
equation to a smaller $x$ a nonvanishing neutralization scale is generated
dynamically. In this sense, the modified MV model should be viewed as
including some of the effects of the low $x$ evolution. Note that we used the
model just for the illustration. In the next section, we will use the MV model
accompanied by the JIMWLK evolution to compute this observable.

The figure 2 demonstrates the presence of the peak at the zero relative angle,
with the peak strength being the strongest for
$|{\underline{k}}_{1}|\approx|{\underline{k}}_{2}|$. We expect that the
numerics of the small x evolution will lead to the transition from the pattern
of the correlation in Fig. 2(a) to the one in Fig. 2(b).

(a) $\alpha_{s}Y=0.0$

(b) $\alpha_{s}Y=0.4$

(c) $\alpha_{s}Y=0.8$

(d) $\alpha_{s}Y=1.0$

Figure 3: Incoherent Diffractive production as a function of relative angle
$k_{1}=7Q^{0}_{s}$, $\epsilon_{f}=2Q^{0}_{s}$.

Note that the Wilson line correlator entering into the inclusive cross-section
to the first nontrivial order is

$\displaystyle\mathcal{N}_{\text{inclusive}}\approx$
$\displaystyle\frac{g^{2}}{N_{c}}\text{Tr}(t^{a}t^{b})\bigg{(}\langle\alpha_{a}({\underline{x}}_{1})\alpha_{b}({\underline{x}}_{1}^{\prime})\rangle-\langle\alpha_{a}({\underline{x}}_{1})\alpha_{b}({\underline{x}}_{2}^{\prime})\rangle-\langle\alpha_{a}({\underline{x}}_{2})\alpha_{b}({\underline{x}}_{1}^{\prime})\rangle+\langle\alpha_{a}({\underline{x}}_{2})\alpha_{b}({\underline{x}}_{2}^{\prime})\rangle\bigg{)}$
(42)

or

$\displaystyle\mathcal{N}_{\text{inclusive}}\approx$
$\displaystyle\frac{C_{f}g^{4}\bar{\mu}^{2}}{2}\bigg{(}L({\underline{x}}_{1}-{\underline{x}}_{1}^{\prime})-L({\underline{x}}_{1}-{\underline{x}}_{2}^{\prime})-L({\underline{x}}_{2}-{\underline{x}}_{1}^{\prime})+L({\underline{x}}_{2}-{\underline{x}}_{2}^{\prime})\bigg{)}\,.$
(43)

As this expression demonstrates the inclusive production has a non-zero
contribution already at order of $\bar{\mu}^{2}$. This contribution is not
sensitive to quantum gluon correlations. Now at $\bar{\mu}^{4}$ order the
correlations including Bose Enhancement in the target are present, however
they do not result to a peak structure in the angular correlation function at
zero angle. The reason why can be understood by considering the diagram in
Fig. 1(c). The diffractive contribution from this diagram is proportional to a
quadratic combination of $L$ times
$\delta_{ab^{\prime}}\delta_{ba^{\prime}}\left[\frac{1}{N_{c}}\operatorname{tr}(t^{a}t^{b})\right]\left[\frac{1}{N_{c}}\operatorname{tr}(t^{a^{\prime}}t^{b^{\prime}})\right]=\frac{C_{f}}{2N_{c}}$,
while for the inclusive production we have the same combination of $L$
functions times
$\delta_{ab^{\prime}}\delta_{ba^{\prime}}\left[\frac{1}{N_{c}}\operatorname{tr}(t^{a}t^{a^{\prime}}t^{b^{\prime}}t^{b})\right]=-\frac{C_{f}}{2N_{c}}$.
We see that the contribution of the same Bose enhancement term but with the
sign flipped due to color algebra. We thus conclude that Bose enhancement,
which gives an enhancement of the zero relative angle peak
(${\underline{k}}_{1}\approx{\underline{k}}_{2}$) in diffractive production,
gives a dip in inclusive! The dip on the background of falling back-to-back
peak is impossible to distinguish in a real experiment and this is why we did
not pursue its numerical calculation in what follows.

##  V Beyond dilute approximation and small x evolution

The analysis performed in the previous section relied upon the dilute
approximation. In this section, we present the results obtained by numerical
simulations of the full MV model with subsequent small-x evolution computed
using leading-order fixed coupling JIMWLK equation.

JIMWLK is the renormalization group evolution towards smaller x, as
conveniently parametrized by $\alpha_{s}Y=\alpha_{s}\ln\frac{1}{x}$. We
initialize the evolution with the MV model, generalized to include the color
neutralization scale $m_{0}=1/4Q^{0}_{s}$ in (40), where $Q_{s}^{0}$ is the
saturation momentum at zero rapidity. We follow the widely accepted definition
of the saturation scale in the numerical calculation through the fundamental
dipole scattering matrix $S(r_{s}=\sqrt{2}/Q_{s})=\exp(-1/2)$.

As discussed above, the conventional MV model does not incorporate color
neutralization. However the JIMWLK evolution generates a color neutralization
scale dynamically even if it is not present in the initial condition. This
scale turns out to be proportional to the saturation momentum with the
proportionality coefficient of order unity, see Ref. [68]. With this in mind
we include such a scale in the initial condition as well albeit with a smaller
proportionality factor.

Numerical procedures for MV and JIMWLK are identical to those described in
Ref. [69] and Ref. [70] correspondingly. As in the previous section, we fixed
the value of $z$ (although it can be chosen arbitrary), in order to compute
the correlations functions defined in Eq. (39). The details of the numerical
calculation of the diffractive observable are presented in Appendix B.

The results are presented in the plots, see Figs. 3 and 4.

Figure 4: Incoherent Diffractive production for collinear configuration
$\Delta\phi=0$ as a function o f $k_{2}$: $k_{1}=7Q^{0}_{s}$,
$\epsilon_{f}=2Q^{0}_{s}$.

##  VI Conclusions

Let us discuss our numerical results.

First consider Fig. 2. The correlator plotted exhibits the features expected
from the Bose enhancement. Both, for $m=0$ and $m=Q_{s}$ the correlator has a
maximum at zero angle as long as the momenta of the two gluons are close to
each other. The maximum is most pronounced when $k_{1}=k_{2}$ and disappears
when the ratio of the two momenta becomes roughly 1.5. Introduction of the
color neutralization scale makes the maximum more robust, both in terms of it
being more pronounced and surviving for wider range of momenta. This feature
was observed also in our earlier work on trijet production [45, 46].
Qualitatively, the origin of this behavior is that finite color neutralization
scale suppresses gluons with small transverse momenta. These gluons mostly
contribute to the correlator at $\Delta\phi=\pi$, and suppressing them makes
the peak at $\Delta\phi=0$ more pronounced. Thus the presence of a large
saturation (or color neutralization) scale is a welcome feature from the
perspective of observing the Bose enhancement effects.

Moving on to Fig. 3 we see the confirmation of the same trend. Evolution in
energy generates color neutralization scale [68]. This scale quickly grows
beyond $Q_{s}^{0}$ present in the initial condition. As a result evolving from
$\alpha_{s}Y=0$ to $\alpha_{s}Y=.4$ significantly increases the Bose
enhancement signal. This is clearly seen moving from Fig. 3(a) to Fig. 3(b).
Interestingly, further evolution to $\alpha_{s}Y=.8$ and beyond changes the
picture very little. Although there is marginal enhancement of the signal, it
is much less noticeable than the initial change starting with the MV model.
This suggests that there is importance not only to the value of the saturation
scale, but also to the shape of the gluon distribution in the hadron.

To summarize, we have shown that one can probe Bose enhancement with a simpler
observable than discussed earlier - the diffractive dijet production. To do
that one needs to go away from the correlation limit into regime where the
total momentum and the momentum imbalance of the two jets are close.

###### Acknowledgements.

A.K. is supported by the NSF Nuclear Theory grant 2208387. This material is
based upon work supported by the U.S. Department of Energy, Office of Science,
Office of Nuclear Physics through the Contract No. DE-SC0020081 (V.S.) and the
Saturated Glue (SURGE) Topical Collaboration. A.K. and V.S. thank the
Binational Science Foundation grant #2021789 for support. V.S. thanks the
ExtreMe Matter Institute for partial support and A. Andronic for hospitality
at Physics Department of Muenster University.

##  Appendix A Useful sums

Using the identity

$\displaystyle(n+m)!=\int_{0}^{\infty}ste^{-t}t^{n+m}$ (44)

we get

$\displaystyle\sum_{n,m}n\rho_{n,m,n,m}$
$\displaystyle=(1-R)\sum_{n,m}n\frac{(n+m)!}{n!m!}\left(\frac{R}{2}\right)^{n+m}=(1-R)\int_{0}^{\infty}dte^{-t}\sum_{n}n\frac{1}{n!}\left(\frac{Rt}{2}\right)^{n}\sum_{m}\frac{1}{m!}\left(\frac{Rt}{2}\right)^{m}$
$\displaystyle=\frac{(1-R)R}{2}\int_{0}^{\infty}dt\,t\,e^{-t(1-R)}=\frac{R}{2(R-1)}=\frac{g^{2}\mu^{2}}{q^{2}}\,.$
(45)

The sum

$\displaystyle\sum_{n,m}\sqrt{m+1}\sqrt{n+1}\rho_{n,m,n+1,m+1}$
$\displaystyle=(1-R)\sum_{n,m}\sqrt{m+1}\sqrt{n+1}\frac{(n+m+1)!}{\sqrt{n!m!(n+1)!(m+1)!}}\left(\frac{R}{2}\right)^{n+m+1}$
$\displaystyle=(1-R)\sum_{n,m}\frac{(n+m+1)!}{n!m!}\left(\frac{R}{2}\right)^{n+m+1}=(1-R)\sum_{n,m}n\frac{(n+m)!}{n!m!}\left(\frac{R}{2}\right)^{n+m}$
(46)

reduces to Eq. (A).

##  Appendix B Numerical evaluation of the incoherent diffractive cross-
section

The main difficulty in computing the expectation value for both diffractive
and inclusive cross section is in computing the integrals of the term
involving four Wilson lines

$\displaystyle T=$
$\displaystyle\int\frac{\text{d}^{2}{\underline{x}}_{1}}{(2\pi)^{2}}\frac{\text{d}^{2}{\underline{x}}_{1}^{\prime}}{(2\pi)^{2}}\frac{\text{d}^{2}{\underline{x}}_{2}}{(2\pi)^{2}}\frac{\text{d}^{2}{\underline{x}}_{2}^{\prime}}{(2\pi)^{2}}e^{-i{\underline{k}}_{1}\cdot({\underline{x}}_{1}-{\underline{x}}_{1}^{\prime})}e^{-i{\underline{k}}_{2}\cdot({\underline{x}}_{2}-{\underline{x}}_{2}^{\prime})}K_{0}(\epsilon_{f}|{\underline{x}}_{1}-{\underline{x}}_{2}|)K_{0}(\epsilon_{f}|{\underline{x}}_{1}^{\prime}-{\underline{x}}_{2}^{\prime}|)$
$\displaystyle\times\left[V^{\dagger}({\underline{x}}_{2})V({\underline{x}}_{1})\right]_{ij}\left[V^{\dagger}({\underline{x}}^{\prime}_{1})V({\underline{x}}^{\prime}_{2})\right]_{kl}\,.$
(47)

To proceed define the amplitude

$\displaystyle M_{ij}({\underline{k}}_{1},{\underline{k}}_{2})=$
$\displaystyle\int\frac{\text{d}^{2}{\underline{x}}_{1}}{(2\pi)^{2}}\frac{\text{d}^{2}{\underline{x}}_{2}}{(2\pi)^{2}}e^{-i{\underline{k}}_{1}\cdot{\underline{x}}_{1}-i{\underline{k}}_{2}{\underline{x}}_{2}}K_{0}(\epsilon_{f}|{\underline{x}}_{1}-{\underline{x}}_{2}|)\left[V^{\dagger}({\underline{x}}_{2})V({\underline{x}}_{1})\right]_{ij}$
(48)

which can be efficiently evaluated using the fast Fourier transform.

Once the amplitude is known, we can find the compination

$\displaystyle
T=M_{ij}({\underline{k}}_{1},{\underline{k}}_{2})M_{kl}(-{\underline{k}}_{2},-{\underline{k}}_{1})$
(49)

and thus for the diffractive we simply

$\displaystyle T={\rm tr}M({\underline{k}}_{1},{\underline{k}}_{2})\,{\rm
tr}M(-{\underline{k}}_{2},-{\underline{k}}_{1})$ (50)

which requires storing a four dimensional matrix ${\rm
tr}M({\underline{k}}_{1},{\underline{k}}_{2})$ as enumerates by the momentum
components. For inclusive one has to store full
$M_{ij}({\underline{k}}_{1},{\underline{k}}_{2})$. That is computing the
inclusive cross section requires a factor of $N_{c}^{2}-1$ more operating
memory and thus is more challenging.

Finally, the diffractive cross section boils down to evaluating the averages
over the target proportional to

$\displaystyle\langle{\rm tr}M({\underline{k}}_{1},{\underline{k}}_{2})\,{\rm
tr}M(-{\underline{k}}_{2},-{\underline{k}}_{1})\rangle-\langle{\rm
tr}M({\underline{k}}_{1},{\underline{k}}_{2})\,\rangle\langle{\rm
tr}M(-{\underline{k}}_{2},-{\underline{k}}_{1})\rangle\,.$ (51)

## References

  * [1] E. C. Aschenauer, S. Fazio, J. H. Lee, H. Mantysaari, B. S. Page, B. Schenke, T. Ullrich, R. Venugopalan, and P. Zurita.  The electron–ion collider: assessing the energy dependence of key measurements.  Rept. Prog. Phys., 82(2):024301, 2019.
  * [2] L. V. Gribov, E. M. Levin, and M. G. Ryskin.  Semihard Processes in QCD.  Phys. Rept., 100:1–150, 1983.
  * [3] Alfred H. Mueller and Jian-wei Qiu.  Gluon Recombination and Shadowing at Small Values of x.  Nucl. Phys. B, 268:427–452, 1986.
  * [4] Alfred H. Mueller.  Parton saturation at small x and in large nuclei.  Nucl. Phys. B, 558:285–303, 1999.
  * [5] Alfred H. Mueller.  Parton saturation: An Overview.  In Cargese Summer School on QCD Perspectives on Hot and Dense Matter, pages 45–72, 11 2001.
  * [6] Edmond Iancu and Raju Venugopalan.  The Color glass condensate and high-energy scattering in QCD, pages 249–3363.  3 2003.
  * [7] Francois Gelis, Edmond Iancu, Jamal Jalilian-Marian, and Raju Venugopalan.  The Color Glass Condensate.  Ann. Rev. Nucl. Part. Sci., 60:463–489, 2010.
  * [8] Yuri V. Kovchegov and Eugene Levin.  Quantum Chromodynamics at High Energy, volume 33.  Oxford University Press, 2013.
  * [9] Heikki Mäntysaari, Niklas Mueller, Farid Salazar, and Björn Schenke.  Multigluon Correlations and Evidence of Saturation from Dijet Measurements at an Electron-Ion Collider.  Phys. Rev. Lett., 124(11):112301, 2020.
  * [10] Cyrille Marquet.  Forward inclusive dijet production and azimuthal correlations in p(A) collisions.  Nucl. Phys. A, 796:41–60, 2007.
  * [11] Daniel Boer, Stanley J. Brodsky, Piet J. Mulders, and Cristian Pisano.  Direct Probes of Linearly Polarized Gluons inside Unpolarized Hadrons.  Phys. Rev. Lett., 106:132001, 2011.
  * [12] Fabio Dominguez, Cyrille Marquet, Bo-Wen Xiao, and Feng Yuan.  Universality of Unintegrated Gluon Distributions at small x.  Phys. Rev. D, 83:105005, 2011.
  * [13] Fabio Dominguez, Jian-Wei Qiu, Bo-Wen Xiao, and Feng Yuan.  On the linearly polarized gluon distributions in the color dipole model.  Phys. Rev. D, 85:045003, 2012.
  * [14] Andreas Metz and Jian Zhou.  Distribution of linearly polarized gluons inside a large nucleus.  Phys. Rev. D, 84:051503, 2011.
  * [15] Cristian Pisano, Daniel Boer, Stanley J. Brodsky, Maarten G. A. Buffing, and Piet J. Mulders.  Linear polarization of gluons and photons in unpolarized collider experiments.  JHEP, 10:024, 2013.
  * [16] A. H. Mueller, Bo-Wen Xiao, and Feng Yuan.  Sudakov double logarithms resummation in hard processes in the small-x saturation formalism.  Phys. Rev. D, 88(11):114010, 2013.
  * [17] P. Kotko, K. Kutak, C. Marquet, E. Petreska, S. Sapeta, and A. van Hameren.  Improved TMD factorization for forward dijet production in dilute-dense hadronic collisions.  JHEP, 09:106, 2015.
  * [18] Adrian Dumitru, Tuomas Lappi, and Vladimir Skokov.  Distribution of Linearly Polarized Gluons and Elliptic Azimuthal Anisotropy in Deep Inelastic Scattering Dijet Production at High Energy.  Phys. Rev. Lett., 115(25):252301, 2015.
  * [19] Adrian Dumitru and Vladimir Skokov.  $cos(4\phi$) azimuthal anisotropy in small-$x$ DIS dijet production beyond the leading power TMD limit.  Phys. Rev. D, 94(1):014030, 2016.
  * [20] A. van Hameren, P. Kotko, K. Kutak, C. Marquet, E. Petreska, and S. Sapeta.  Forward di-jet production in p+Pb collisions in the small-x improved TMD factorization framework.  JHEP, 12:034, 2016.  [Erratum: JHEP 02, 158 (2019)].
  * [21] C. Marquet, E. Petreska, and C. Roiesnel.  Transverse-momentum-dependent gluon distributions from JIMWLK evolution.  JHEP, 10:065, 2016.
  * [22] P. Kotko, K. Kutak, S. Sapeta, A. M. Stasto, and M. Strikman.  Estimating nonlinear effects in forward dijet production in ultra-peripheral heavy ion collisions at the LHC.  Eur. Phys. J. C, 77(5):353, 2017.
  * [23] Cyrille Marquet, Claude Roiesnel, and Pieter Taels.  Linearly polarized small-$x$ gluons in forward heavy-quark pair production.  Phys. Rev. D, 97(1):014004, 2018.
  * [24] Adrian Dumitru, Vladimir Skokov, and Thomas Ullrich.  Measuring the Weizsäcker-Williams distribution of linearly polarized gluons at an electron-ion collider through dijet azimuthal asymmetries.  Phys. Rev. C, 99(1):015204, 2019.
  * [25] Edmond Iancu and Yair Mulian.  Forward trijet production in proton–nucleus collisions.  Nucl. Phys. A, 985:66–127, 2019.
  * [26] Spencer R. Klein and Heikki Mäntysaari.  Imaging the nucleus with high-energy photons.  Nature Rev. Phys., 1(11):662–674, 2019.
  * [27] Edmond Iancu and Yair Mulian.  Forward dijets in proton-nucleus collisions at next-to-leading order: the real corrections.  JHEP, 03:005, 2021.
  * [28] Yoshitaka Hatta, Bo-Wen Xiao, Feng Yuan, and Jian Zhou.  Azimuthal angular asymmetry of soft gluon radiation in jet production.  Phys. Rev. D, 104(5):054037, 2021.
  * [29] Renaud Boussarie, Heikki Mäntysaari, Farid Salazar, and Björn Schenke.  The importance of kinematic twists and genuine saturation effects in dijet production at the Electron-Ion Collider.  JHEP, 09:178, 2021.
  * [30] Hirotsugu Fujii, Cyrille Marquet, and Kazuhiro Watanabe.  Comparison of improved TMD and CGC frameworks in forward quark dijet production.  JHEP, 12:181, 2020.
  * [31] Benjamin Rodriguez-Aguilar, D. N. Triantafyllopoulos, and S. Y. Wei.  Incoherent diffractive dijet production in electron DIS off nuclei at high energy.  Phys. Rev. D, 107(11):114007, 2023.
  * [32] Paul Caucal, Farid Salazar, and Raju Venugopalan.  Dijet impact factor in DIS at next-to-leading order in the Color Glass Condensate.  JHEP, 11:222, 2021.
  * [33] Paul Caucal, Farid Salazar, Björn Schenke, Tomasz Stebel, and Raju Venugopalan.  Back-to-back inclusive dijets in DIS at small x: gluon Weizsäcker-Williams distribution at NLO.  JHEP, 08:062, 2023.
  * [34] Paul Caucal, Farid Salazar, Björn Schenke, Tomasz Stebel, and Raju Venugopalan.  Back-to-back inclusive dijets in DIS at small $x$: Complete NLO results and predictions.  7 2023.
  * [35] Jochen Bartels, H. Jung, and M. Wusthoff.  Quark - anti-quark gluon jets in DIS diffractive dissociation.  Eur. Phys. J. C, 11:111–125, 1999.
  * [36] Tolga Altinoluk, Néstor Armesto, Guillaume Beuf, and Amir H. Rezaeian.  Diffractive Dijet Production in Deep Inelastic Scattering and Photon-Hadron Collisions in the Color Glass Condensate.  Phys. Lett. B, 758:373–383, 2016.
  * [37] Yoshitaka Hatta, Bo-Wen Xiao, and Feng Yuan.  Probing the Small- x Gluon Tomography in Correlated Hard Diffractive Dijet Production in Deep Inelastic Scattering.  Phys. Rev. Lett., 116(20):202301, 2016.
  * [38] Yoshikazu Hagiwara, Yoshitaka Hatta, Roman Pasechnik, Marek Tasevsky, and Oleg Teryaev.  Accessing the gluon Wigner distribution in ultraperipheral $pA$ collisions.  Phys. Rev. D, 96(3):034009, 2017.
  * [39] Heikki Mäntysaari, Niklas Mueller, and Björn Schenke.  Diffractive Dijet Production and Wigner Distributions from the Color Glass Condensate.  Phys. Rev. D, 99(7):074004, 2019.
  * [40] Farid Salazar and Björn Schenke.  Diffractive dijet production in impact parameter dependent saturation models.  Phys. Rev. D, 100(3):034007, 2019.
  * [41] E. Iancu, A. H. Mueller, and D. N. Triantafyllopoulos.  Probing Parton Saturation and the Gluon Dipole via Diffractive Jet Production at the Electron-Ion Collider.  Phys. Rev. Lett., 128(20):202001, 2022.
  * [42] Yoshitaka Hatta, Bo-Wen Xiao, and Feng Yuan.  Semi-inclusive diffractive deep inelastic scattering at small x.  Phys. Rev. D, 106(9):094015, 2022.
  * [43] G. Beuf, H. Hänninen, T. Lappi, Y. Mulian, and H. Mäntysaari.  Diffractive deep inelastic scattering at NLO in the dipole picture: The qq¯g contribution.  Phys. Rev. D, 106(9):094014, 2022.
  * [44] E. Iancu, A. H. Mueller, D. N. Triantafyllopoulos, and S. Y. Wei.  Gluon dipole factorisation for diffractive dijets.  JHEP, 10:103, 2022.
  * [45] Alex Kovner, Ming Li, and Vladimir V. Skokov.  Probing Gluon Bose Correlations in Deep Inelastic Scattering.  Phys. Rev. Lett., 128(18):182003, 2022.
  * [46] Alex Kovner, Ming Li, and Vladimir V. Skokov.  Probing gluon Bose correlations in nuclear wave function in deep inelastic scattering.  Phys. Rev. D, 107(11):114032, 2023.
  * [47] Adrian Dumitru, Kevin Dusling, Francois Gelis, Jamal Jalilian-Marian, Tuomas Lappi, and Raju Venugopalan.  The Ridge in proton-proton collisions at the LHC.  Phys. Lett. B, 697:21–25, 2011.
  * [48] Tolga Altinoluk, Néstor Armesto, Guillaume Beuf, Alex Kovner, and Michael Lublinsky.  Bose enhancement and the ridge.  Phys. Lett. B, 751:448–452, 2015.
  * [49] Alex Kovner and Vladimir V. Skokov.  Bose enhancement, the Liouville effective action and the high multiplicity tail in p-A collisions.  Phys. Rev. D, 98(1):014004, 2018.
  * [50] Jonathan R. Gaunt and W. James Stirling.  Double Parton Distributions Incorporating Perturbative QCD Evolution and Momentum and Quark Number Sum Rules.  JHEP, 03:005, 2010.
  * [51] B. Blok, Yu. Dokshitzer, L. Frankfurt, and M. Strikman.  The Four jet production at LHC and Tevatron in QCD.  Phys. Rev. D, 83:071501, 2011.
  * [52] Markus Diehl, Daniel Ostermeier, and Andreas Schafer.  Elements of a theory for multiparton interactions in QCD.  JHEP, 03:089, 2012.  [Erratum: JHEP 03, 001 (2016)].
  * [53] B. Blok, Yu. Dokshitzer, L. Frankfurt, and M. Strikman.  Perturbative QCD correlations in multi-parton collisions.  Eur. Phys. J. C, 74:2926, 2014.
  * [54] Markus Diehl and Jonathan R. Gaunt.  Double parton scattering theory overview.  Adv. Ser. Direct. High Energy Phys., 29:7–28, 2018.
  * [55] V. V. Sudakov.  Vertex parts at very high-energies in quantum electrodynamics.  Sov. Phys. JETP, 3:65–71, 1956.
  * [56] Larry D. McLerran and Raju Venugopalan.  Gluon distribution functions for very large nuclei at small transverse momentum.  Phys. Rev. D, 49:3352–3355, 1994.
  * [57] Larry D. McLerran and Raju Venugopalan.  Computing quark and gluon distribution functions for very large nuclei.  Phys. Rev. D, 49:2233–2241, 1994.
  * [58] Tolga Altinoluk, Guillaume Beuf, Alina Czajka, and Arantxa Tymowska.  DIS dijet production at next-to-eikonal accuracy in the CGC.  Phys. Rev. D, 107(7):074016, 2023.
  * [59] Giovanni Antonio Chirilli.  Sub-eikonal corrections to scattering amplitudes at high energy.  JHEP, 01:118, 2019.
  * [60] Florian Cougoulic, Yuri V. Kovchegov, Andrey Tarasov, and Yossathorn Tawabutr.  Quark and gluon helicity evolution at small x: revised and updated.  JHEP, 07:095, 2022.
  * [61] Ming Li.  Small x physics beyond eikonal approximation: an effective Hamiltonian approach.  JHEP, 07:158, 2023.
  * [62] Alex Kovner and Michael Lublinsky.  Entanglement entropy and entropy production in the Color Glass Condensate framework.  Phys. Rev. D, 92(3):034016, 2015.
  * [63] Haowu Duan, Candost Akkaya, Alex Kovner, and Vladimir V. Skokov.  Entanglement, partial set of measurements, and diagonality of the density matrix in the parton model.  Phys. Rev. D, 101(3):036017, 2020.
  * [64] Haowu Duan, Alex Kovner, and Vladimir V. Skokov.  Gluon quasiparticles and the CGC density matrix.  Phys. Rev. D, 105(5):056009, 2022.
  * [65] T. Lappi, B. Schenke, S. Schlichting, and R. Venugopalan.  Tracing the origin of azimuthal gluon correlations in the color glass condensate.  JHEP, 01:061, 2016.
  * [66] Edmond Iancu, Kazunori Itakura, and Larry McLerran.  A Gaussian effective theory for gluon saturation.  Nucl. Phys. A, 724:181–222, 2003.
  * [67] Larry McLerran and Vladimir Skokov.  Finite Numbers of Sources, Particle Correlations and the Color Glass Condensate.  Nucl. Phys. A, 947:142–154, 2016.
  * [68] Haowu Duan, Vladimir Skokov, and Cameron Stephens.  Color neutralization scale at small x.  2023\.
  * [69] T. Lappi.  Wilson line correlator in the MV model: Relating the glasma to deep inelastic scattering.  Eur. Phys. J. C, 55:285–292, 2008.
  * [70] T. Lappi and H. Mäntysaari.  On the running coupling in the JIMWLK equation.  Eur. Phys. J. C, 73(2):2307, 2013.