# Richardson–Lucy deconvolution with a spatially Variant point-spread function
of Chandra: Supernova Remnant Cassiopeia A as an Example

Yusuke Sakai Department of Physics, Rikkyo University, Toshima-Ku, Tokyo,
171-8501, Japan  Shinya Yamada Department of Physics, Rikkyo University,
Toshima-Ku, Tokyo, 171-8501, Japan  Toshiki Sato Department of Physics, Rikkyo
University, Toshima-Ku, Tokyo, 171-8501, Japan  Department of Physics, School
of Science and Technology, Meiji University, 1-1-1 Higashi Mita, Tama-ku,
Kawasaki, Kanagawa 214-8571, Japan  Ryota Hayakawa Department of Physics,
Rikkyo University, Toshima-Ku, Tokyo, 171-8501, Japan  International Center
for Quantum-field Measurement Systems for Studies of the Universe and
Particles (QUP), KEK, 1-1 Oho, Tsukuba, Ibaraki 305-0801, Japan  Ryota
Higurashi Department of Physics, Rikkyo University, Toshima-Ku, Tokyo,
171-8501, Japan  Nao Kominato Department of Physics, Rikkyo University,
Toshima-Ku, Tokyo, 171-8501, Japan

###### Abstract

Richardson–Lucy (RL) deconvolution is one of the classical methods widely used
in X-ray astronomy and other areas. Amid recent progress in image processing,
RL deconvolution still leaves much room for improvement under a realistic
situations. One direction is to include the positional dependence of a point-
spread function (PSF), so-called RL deconvolution with a spatially variant PSF
(RLsv). Another is the method of estimating a reliable number of iterations
and their associated uncertainties. We developed a practical method that
incorporates the RLsv algorithm and the estimation of uncertainties. As a
typical example of bright and high-resolution images, the Chandra X-ray image
of the supernova remnant Cassiopeia A was used in this paper. RLsv
deconvolution enables us to uncover the smeared features in the
forward/backward shocks and jet-like structures. We constructed a method to
predict the appropriate number of iterations by using statistical fluctuation
of the observed images. Furthermore, the uncertainties were estimated by error
propagation from the last iteration, which was phenomenologically tested with
the observed data. Thus, our method is a practically efficient framework to
evaluate the time evolution of the remnants and their fine structures embedded
in high-resolution X-ray images.

Astronomy data analysis (1858), Astronomy image processing (2306), High
angular resolution (2167), X-ray astronomy (1810)

##  1 Introduction

Imaging analysis is critically important for studying diffuse celestial
sources. X-ray astronomy, starting with the first space application of X-ray
CCD in ASCA (Burke et al., 1994), has delivered detailed images of various
celestial objects; e.g., supernova remnants such as SN1006 (Bamba et al.,
2003) and Cassiopeia A (hereafter Cas A; Hwang et al., 2004), and galaxy
clusters such as A2142 (Markevitch et al., 2000). Since the X-ray mirrors used
in Chandra are the largest and most precisely built, exceeding the angular
resolution of Chandra is considered to be challenging. Therefore, enhancing
the technique of imaging analysis has been an essential direction to utilize
the highest spatial resolution and data accumulated over decades.

X-rays are collected primarily by total reflection from the surface of an
X-ray mirror, therefore the response function for the distribution of focused
X-rays, called the point-spread function (PSF), is nearly energy independent.
On the condition that a PSF is independent of incoming photon energy and the
position of the focal plane, a reverse calculation of a convolution of PSF,
so-called image deconvolution (see the review on deconvolution in astronomy by
Starck et al. (2002)), is highly simplified. There are various deconvolution
methods proposed by assuming that a PSF is constant during the deconvolution
process, such as the deconvolution of Suzaku XIS (Sugizaki et al., 2009). One
of the latest examples is the image restoration algorithm Expectation via
Markov chain Monte Carlo (Esch et al., 2004). It is applied to the double
active galactic nuclei in NGC 6240 (e.g., Fabbiano et al., 2020; Paggi et al.,
2022), succeeding in finely resolving the two cores. Similarly, a classical
method, Richardson–Lucy (RL) deconvolution proposed by Richardson (1972) and
Lucy (1974), is often used (e.g., Grefenstette et al., 2015; Thimmappa et al.,
2020; Sobolenko et al., 2022).

The choice of method depends on the trade-off between accuracy and
computational cost. Relaxing the condition that a PSF is positional-
independent and/or energy-independent, the deconvolution methods increase the
complexity of the calculation. RL deconvolution is one of the simplified
methods but still has room for improvement in practical situations. In gamma-
ray astronomy, a PSF can change by one order of magnitude with energy and
incident angle; it is calculated for each event, e.g., RL algorithm optimized
for Fermi-LAT and EGRET (Tajima et al., 2007). In contrast, as the number of
photons is much larger in X-ray astronomy, event-by-event reconstruction is
less practical; image-based reconstruction thus can be the first choice in
X-rays. However, there are few studies on extending the RL method, especially
their application to diffuse sources obtained by Chandra. We therefore
explored its applicability to the Chandra data and considered the associated
systematic errors.

In this paper, we implement RL deconvolution with a spatially variant PSF
(RLsv) algorithm, assuming it to be used for Chandra images. Section 2
describes the principle of the RLsv method. One of the technical difficulties
is reducing computational cost in calculating PSFs. This is solved by
decimating the sampling interval of PSFs, while a side-effect is discussed in
Section 5.1. Section 3 presents an example of its application to a diffuse
source observed by Chandra. We apply the RLsv method to the supernova remnant
of Cas A as an example, because Cas A is bright and extended over the entire
field of view of the ACIS detector, which would be the best target for the
first application. The remnant is intensively studied because of its unique
structure and evolution, e.g., the velocities and thickness of shocked
filaments (Patnaude & Fesen, 2009; Sato et al., 2018; Tsuchioka et al., 2022),
where the method can contribute to advancing our understanding of the
phenomena. In Section 4, we propose a reliable number of stop iterations and
uncertainties of the method. We develop the method to estimate the number of
convergent iterations by generating fluctuations due to statistical errors
during iteration. Furthermore, the uncertainty on the RLsv-deconvolved image
is estimated by using the law of error propagation (e.g., Ku, 1966). As a
result, filaments and ambiguous structures of Cas A are deconvolved to be
sharper with some knowledge of the statistical uncertainties.

Table 1: Basic information on the Chandra Observations of Cas A Used in this Paper Obs. ID | Obs. Start | Exp. Time | detector | R.A. | Decl. | Roll
---|---|---|---|---|---|---
| yyyy mmm dd | (ks) |  | (deg) | (deg) | (deg)
4636 | 2004 Apr 20 | 143.48 | ACIS-S | 350.9129 | 58.8412 | 49.7698
4637 | 2004 Apr 22 | 163.50 | ACIS-S | 350.9131 | 58.8414 | 49.7665
4639 | 2004 Apr 25 | 79.05 | ACIS-S | 350.9132 | 58.8415 | 49.7666
5319 | 2004 Apr 18 | 42.26 | ACIS-S | 350.9127 | 58.8411 | 49.7698
5196 | 2004 Feb 8 | 49.53 | ACIS-S | 350.9129 | 58.7933 | 325.5035

##  2 Method

###  2.1 RL Deconvolution

The RL algorithm iteratively estimates a true image from an observed image
using Bayesian inference. It generally assumes that the PSF does not change
with a position in the image. The RL algorithm is expressed by

$W_{i}^{(r+1)}=W_{i}^{(r)}\sum_{k}\frac{P_{ik}H_{k}}{\sum_{j}P_{jk}W_{j}^{(r)}},$
(1)

where $i$ and $j$ are mapping the image in the sky, and $k$ is mapping the
image on the detector. The indices of the summation run through all the
pixels. $W^{(r)}$ is the restored image after $r$ iterations, and $H$ is the
observed image on the ACIS detector. $P_{jk}$ is the probability that a photon
emitted in sky $W$ bin $j$ is measured in data space $H$ bin $k$, or
$P(H_{k}|W_{j})$.

###  2.2 RL with a Spatially Variant PSF

Previous Chandra image deconvolution approaches (e.g., Thimmappa et al., 2020;
Sobolenko et al., 2022) used a simplified approximation for the $P_{jk}$
values, i.e., they used the same PSF for each $j$ bin. Here we assume that the
PSF changes as a function of the off-axis angle and the roll angle. As a
consequence, the Chandra RL algorithm is extended. The formula for RLsv is
obtained by rewriting Equation (1) as

$W_{i}^{(r+1)}=W_{i}^{(r)}\sum_{k}\frac{P_{iik}H_{k}}{\sum_{j}P_{jjk}W_{j}^{(r)}}.$
(2)

$P_{jjk}$ refers to a PSF at a position of $j$ (first index) which returns a
probability that an event emitted at $W_{j}$ (second index) is observed at
$H_{k}$ (third index), or $P_{j}(H_{k}|W_{j})$. Computational cost and memory
requirements need to be minimized for calculating the third-order tensor of
the PSF, which is a distinctive feature of the RLsv algorithm. When $H$ is
corrected for slight differences among pixels in effective area and exposures,
its normalization can be chosen arbitrarily. Here we use $H_{k}=N_{k}/A_{k}$,
where $N_{k}$ is the detector count image, and $A_{k}$ is the Hadamard product
of effective area and exposure time.

###  2.3 RLsv Deconvolution with Total Variation Regularization

There are regularization techniques to enhance the RL method, which are also
readily available for the RLsv method. Among these, total variation (TV)
regularization (Rudin et al., 1992) is effective in handling statistical
errors, which is used in the RL method (Dey et al., 2006). The formula for
RLsv with the regularization is obtained by rewriting Equation (2) as

$W_{i}^{(r+1)}=\frac{W_{i}^{(r)}}{1-\lambda_{\rm{TV}}\textrm{div}\left(\frac{\nabla
W_{i}^{(r)}}{|\nabla
W_{i}^{(r)}|}\right)}\sum_{k}\frac{P_{iik}H_{k}}{\sum_{j}P_{jjk}W_{j}^{(r)}}.$
(3)

The only difference from the RLsv algorithm of Equation (2) is the
regularization term of $1-\lambda_{\rm{TV}}\textrm{div}(\nabla
W_{i}^{(r)}/|\nabla W_{i}^{(r)}|)$, where $\lambda_{\rm{TV}}$ is the
regularization parameter, $\textrm{div}(\cdot)$ is the divergence, and $\nabla
W_{i}^{(r)}$ is the gradient of $W_{i}^{(r)}$. In this paper, we utilize the
parameter of $\lambda_{\rm{TV}}=0.002$, as proposed by Dey et al. (2006).

###  2.4 Comparison to other methods

Deconvolution methods require an understanding of their applicability to a
practical condition, as well as optimization of computation cost and accuracy
(for features of various methods see Naik & Sahu (2013)). The RL method is
well studied and has been used to incorporate regularization (e.g., van Kempen
& van Vliet, 2000; Dey et al., 2006; Yuan et al., 2008; Yongpan et al., 2010)
and the recent trend of deep learning (e.g., Agarwal et al., 2020). For the
Chandra users, RL deconvolution with a single PSF is frequently used because
the method is already implemented as arestore in the Chandra Interactive
Analysis of Observations (CIAO; Fruscione et al., 2006), Chandra’s standard
data processing package.

Figure 1: Cas A image (Obs. ID=4636) and the two-dimensional probabilities of
the point-spread functions (PSFs). The integral of each PSF is normalized to
be 1. The PSF color scale is a fixed range. The location of the optical axis
is indicated with a green cross.

Compared to other methods, the RL method forces the deconvolved image of each
iteration to be non-negative, and its integral value is conserved.
Additionally, the method converges to the maximum likelihood solution for a
Poisson noise distribution (Shepp & Vardi, 1982), which is suitable for
Chandra images with noise from counting statistics. Depending on the
application, it is less prone to ringing artifacts than inverse PSF-based
methods (e.g., Sekko et al., 1999; Neelamani et al., 2004); see the results of
the comparison by Dalitz et al. (2015). According to White (1994), it is
robust against small errors in the PSF.

Figure 2: (a) X-ray image in the 0.5–7.0 keV band of Cas A obtained with
Chandra. (a-1, 2, 3) Enlarged images specified by the colored frames in (a).
(b) Same as (a), but for the RLsv-deconvolved results. The unit of flux in the
images is $\rm{photons~{}cm^{-2}~{}s^{-1}}$.

##  3 Application to Observed Data

###  3.1 Data selection

Because Cas A is a bright and diffuse X-ray source with a moderately large
apparent diameter, it is an ideal target to demonstrate the RLsv method. It
has been observed by Chandra almost every year since 1999. The Chandra data of
Cas A used in this paper are listed in Table 1: ACIS-S observation of 2004 in
Obs. ID=4636, 4637, 4639, and 5319. The image size is $1489\times 1488$
pixels, or $743^{\prime\prime}\times 742^{\prime\prime}$ given a unit pixel of
$0.^{\prime\prime}492$. Data processing and analysis were performed using CIAO
version 4.13. The data were reprocessed from the level 1 event files by
chandra_repro. Since the roll angle and optical axis of the four observations
are almost the same (the maximum difference of the optical axis location is
about 4 unit pixels), all the events were merged into one by merge_obs. The
total exposure time was 428.29 ks.

###  3.2 Generating the PSF of Chandra

The Chandra telescope system consists of four pairs of nested reflecting
surfaces, configured in the Wolter type I geometry. The high energy response
is achieved by coating the mirrors with iridium. It has attained the highest
angular resolution of $0.^{\prime\prime}492$ among existing X-ray telescopes.
Its mirror of Chandra has been extensively calibrated on the ground and in
orbit (Jerius et al., 2000). The Chandra PSF is positional-dependent, mainly
due to aberrations. The RLsv method includes the position dependence of the
PSF, which is useful for highly extended X-ray sources.

Because creating a PSF for each position is computationally expensive, it is
decimated at some intervals. For reference, creating a PSF takes several
seconds, depending on the computational environment and desired accuracy. The
sampling interval of PSFs was chosen to be $35\times 35$ pixels (total of
$43\times 43=1849$). The interval was determined empirically by trying several
different ranges. In general, the PSFs simulated from each observation should
be merged for a precise calculation. Here, the PSF of the Obs. ID of 4636 was
used as a representative since its sampling of the PSF is decimated. The PSFs
at the lattice points were generated by CIAO’s simulate_psf using the Model of
AXAF Response to X-rays (Wise et al., 1997; Davis et al., 2012) at a
monochromatic energy of 2.3 keV. They were applied to the observed image with
energies from 0.5 to 7.0 keV.

Figure 1 shows all the PSFs sampled every $35\times 35$ pixels. The optical
axis is located at the northeast in the image, where the spread of the PSF is
minimum. As the position is away from the optical axis, the tail of the PSF
increases with its gradual shift of the elliptical axis. Although it is a
trade-off with photon statistics, it is effective to run the RLsv method with
the optimal monoenergetic PSF for each of the multiple energy decompositions
(see Figure 6). This is because, at shorter wavelengths, the effect of diffuse
reflection due to the roughness of the mirror surface is not negligible
(Jerius et al., 2000).

###  3.3 Results of the RLsv method

Figure 2(a) is an observed $\sim$400 ks image using the energy range from 0.5
to 7.0 keV, as explained in Section 3.1. We applied the RLsv method to the
image with a sampling interval of each PSF of $35\times 35$ pixels. The number
of iterations is 200. Note that the choice of the iteration number is
discussed in Section 4.1. The result of the RLsv method is presented in Figure
2(b). The unit of flux and its range in Figure 2(b) is the same as in Figure
2(a). The overall structures in the RLsv image become more vivid than the
original ones. The images around the off-axis are significantly improved
compared to those around the optical axis.

To make the differences more precise, we present magnified images of the
original image in Figures 2(a-1, 2, 3) and the RLsv-image in Figures 2(b-1, 2,
3). The three regions represent a sharp filament in the northeast, complicated
filaments in the north, and a slightly diffuse area in the south. The
filamentary structures in the northeast and north become sharper in the RLsv-
image. We will quantify the filament width in detail and discuss the
systematic uncertainties associated with the method in Section 4.2.

##  4 Uncertainty Estimation

###  4.1 Assessment of the reasonable number of iterations

We considered a way of assessing an appropriate number of iterations, which is
one of the issues with the RL method. This is because the method has the
property of excessive amplification of noise as the number of iterations
increases. We propose a method to suppress convergence using statistical
errors during the iteration. The formula of the method is written as

$W_{i}^{(r+1)}=W_{i}^{(r)}\sum_{k}\frac{P_{iik}G(N_{k})/A_{k}}{\sum_{j}P_{jjk}W_{j}^{(r)}}.$
(4)

The only difference from the RLsv algorithm, Equation (2), is the
$G(N_{k})/A_{k}$ term. $N$ ($\rm{counts}$) is the map of detector counts. $A$
($\rm{photons~{}cm^{-2}~{}s^{-1}}$) is the Hadamard product of effective area
and exposure time. $G(N_{k})$ is a random number generator following a Poisson
distribution with a count in the $k$th pixel of $N_{k}$; i.e.,
$G(N_{k})/A_{k}$ is a flux in units of $\rm{photons~{}cm^{-2}~{}s^{-1}}$. The
reason for normalizing $G(N)$ by dividing it with $A$ is to account for the
slight variations in effective area and exposure time among pixels.

The performance of the RLsv algorithm using Equation (4) is compared to that
using Equation (2). The convergence is evaluated by using the mean squared
error (MSE) of the two images: one step before and after iterations. Figure 3
shows the history of the MSE during the iteration. The curve obtained by the
RLsv algorithm, Equation (4) saturates at a certain level, while the other
continues to decrease. The saturation level is caused by the injection of
Poisson fluctuation at each step, which is considered as an indicator of
stopping. In Figure 3, the iteration number of $\sim$30 seems appropriate.

Figure 3: Residuals of the two images before and after iterations for the
entire region of Cas A vs. the number of iterations. The results of RLsv
method with and without statistical errors are plotted as a blue and an orange
line, respectively.

###  4.2 Assessment of image blurredness

We then designed a simplified method for evaluating a certain amount of
confidence. The RLsv-deconvolved image should have a similar amount of
fluctuation accompanying the observation image. The principle of the method is
to propagate errors of the converged RLsv image into the next step. Our choice
of using the last step for the error propagation is just to simplify the task.
Here, each error in the observed image is considered statistically
independent. Assuming only uncertainties on $H_{k}$, using the law of error
propagation, the image uncertainty can be expressed as

$\begin{split}\sigma_{W^{\prime}_{i}}&=\sqrt{\sum_{k}\left[\frac{\partial}{\partial
H_{k}}\left(W_{i}\sum_{k}\frac{P_{iik}H_{k}}{\sum_{j}P_{jjk}W_{j}}\right)\sigma_{H_{k}}\right]^{2}}\\\
&=W_{i}\sqrt{\sum_{k}\left(\frac{P_{iik}}{\sum_{j}P_{jjk}W_{j}}\frac{\sqrt{N_{k}}}{A_{k}}\right)^{2}},\end{split}$
(5)

where $W^{\prime}$ is the image of the next iteration number of any estimated
true image of $W$, and $\sqrt{N_{k}}$ is the statistical error of $N_{k}$.

Figure 4: Images and radial profiles of the southeastern filament of Cas A.
(a) Off-axis image of Obs. ID=4636, 4637, 4639, and 5319. (b) Result of RLsv-
image of (a). (c) On-axis image of Obs. ID=5196. (d) Results of the radial
profile in (a), (b), and (c), radially projected from he central compact
object (CCO) using the fan-shaped regions in green. The horizontal axis
represents the distance from the CCO.

We compared the off-axis RL image with the error of Equation (5) to an on-axis
observation in Figure 4. Figures 4(a) and (b) are the off-axis southeastern
images from Figures 2(a) and (b), respectively. Figure 4(c) shows a
southeastern on-axis image of Obs. ID=5196. The exposure time for the on-axis
observation is $\sim$50 ks, resulting in a larger statistical error compared
to the off-axis observation of $\sim$400 ks. The fan-shaped regions in Figures
4(a)–(c), along the filament, are chosen to create the radial profiles, which
is the one-dimensional profile of the photons in each region extending from
the central compact object (CCO) of Cas A toward the outer regions. These
radial profiles were created using dmextract in CIAO. In Figure 4(d), the
framed regions from Figures 4(a)–(c) are color-coded as blue, orange, and
black, respectively. The error bars of the radial profiles in Figure 4(d)
correspond to statistical errors, represented by blue and black, and the
result obtained by applying Equation (5) to the 199th iteration of the RLsv
image, is indicated by orange. From Figure 4(d), the profile of the off-axis
RLsv image agreed with that of the on-axis image within the statistical
errors. This method gives a guideline for a certain level of confidence
associated with the RLsv method.

##  5 Discussion

###  5.1 Enhancement Technique and Possibilities

In this section, further enhancements to the RLsv method are discussed. The
first is to reduce the loss of down-sampling PSFs. The positional dependence
of the PSF does not contain high-frequency components, so the decimation of
PSF sampling should work to some extent. For a small image such as a core plus
jet structure in an active galactic nucleus, keeping a high sampling rate of
the PSFs might be possible. However, for a largely extended source such as a
supernova remnant or galaxy cluster, to minimize the sampling rate is
critically important for practical use. The higher the decimation, the more
emphasized the boundary of the segment. Taking Cas A as an example, the edges
of specific segments clearly appear when the sampling interval is $35\times
35$ pixels. To smooth out the edges, we propose that the PSFs’ boundaries be
randomly selected from nearby PSFs (see more details in the Appendix).

Figure 5: Comparison of the results of the RLsv method without and with total
variation regularization, shown in (a) and (b) respectively.

Second, this method can be developed by incorporating several regularization
methods. We implemented an RLsv method incorporating the TV regularization
expressed in Equation (3). Finally, the RLsv method is naturally applied to
color images. By decomposing observed images into several colors (or energy
bands) and generating PSFs for an appropriate energy in each band, an energy-
dependent RLsv method can be realized.

Figure 6: (a) X-ray RGB (red: 0.2–1.2 keV, green: 1.2–2.0 keV and blue:
2.0–7.0 keV) band images of Cas A obtained with Chandra. (a-1, -2): Enlarged
images specified by the colored frames in (a). (b) Same as (a) except for
RLsv-deconvolved in each energy band. The unit of flux in the images is
$\rm{photons~{}cm^{-2}~{}s^{-1}}$.

We compare the RLsv method with and without the TV regularization. We use the
same image as in Section 3.1. The PSF sampling is $35\times 35$ pixels and the
number of iterations is 200. The PSFs’ boundaries are randomly selected
following the Appendix. Figure 5 presents the enlarged eastern image after
applying the RLsv method to the entire region. Figures 5(a) and (b) show the
results of the RLsv method of Equation (2) and the regularization version of
Equation (3), respectively. The TV regularization preserves the sharp
structure to remain and smoothes out statistical errors. In this way,
regularization can be added to the RLsv method.

We implement the RLsv method including these enhancements and adapt it to the
Cas A observational data described in Section 3.1. Figure 6(a) is the observed
image of Section 3.3 divided into three energies in RGB: 0.2–1.2 keV (red),
1.2–2.0 keV (green), and 2.0–7.0 keV (blue). Cas A is dominated by the thermal
radiation in $\leqslant$4 keV and the nonthermal radiation in $\geqslant$4
keV. We applied the RLsv method with the TV regularization Equation (3) to
each energy image of Figure 6(a) using the appropriate energy of the PSF (red
is 0.92 keV, green is 1.56 keV, and blue is 3.8 keV) based on the official
CIAO page.111https://cxc.cfa.harvard.edu/ciao/why/monochromatic_energy.html
The sampling interval of PSFs is $35\times 35$ pixels. PSF is randomly
selected at the sampling boundaries according to the Appendix. The number of
iterations is 30, according to Section 4.2. The result of the RLsv method is
presented in Figure 6(b). The energy dependence in Figures 6(b-1) and (b-2)
are clearly visible by this method.

###  5.2 Constraint on the Uncertainty

We propose two complementary ways to obtain a guideline on the stop condition
of the RLsv method. One is to obtain a minimum of residuals by inserting
statistical uncertainties into each update during the iteration process. This
gives a rough estimate of the limit of the iterations. The other is to include
the statistical uncertainties in the last step of the iteration. By combining
the two methods, it is possible to derive uncertainties using the errors
obtained by the latter method at the optimal iteration number estimated by the
former. This is a quick and convenient way to derive the systematic
uncertainties associated with the RLsv method.

The systematic uncertainty in the former method is a way of defining the
optimal number of iterations. One easy way is to use the same level of
residuals as shown in Section 4.1. It is intrinsically difficult to
distinguish the signal of the celestial objects from the statistical noise.
This difficulty needs to be overcome by comparing the deconvolved images
without errors and the error-estimated images around the optimal iteration
number recommended by the former method. This method is based on a compromise
between the computational cost and the simplicity of use while keeping a
reasonable statistical error.

##  6 Conclusion

We have improved the processing capability of RL deconvolution by
incorporating the positional dependency of the Chandra PSF. The RLsv method is
applied to the entire region of Cas A with an estimation of its limit and
errors, which are based on the phenomenological method for evaluating a
reasonable number of iterations and uncertainties. It shows that the features
of shock waves and jets are sharper than those measured in the original image,
with a certain amount of knowledge of the associated errors. The RLsv-
deconvolved profile of the off-axis image at the southeastern filament became
shaper and agreed with that of the on-axis observation within the statistical
errors. This method is useful for a detailed diagnosis of other extended X-ray
sources obtained by Chandra.

The code used in this paper is available at doi:10.5281/zenodo.8020557.

We would like to thank the anonymous referee for helpful comments and feedback
on this paper. This research has made use of data obtained from the Chandra
Data Archive and the Chandra Source Catalog, and software provided by the
Chandra X-ray Center (CXC) in the application packages CIAO. This work was
supported by JSPS KAKENHI grant Nos. 20K20527, 22H01272, and 20H01941.

## Appendix
Boundaries of the PSF

Figure 7: Illustration of the $9\times 9$ pixels around the intersection of
the PSF switchover, overlaid with the probability weights on selecting PSFs.
The quadrants are named A, B, C, and D for illustrative purposes. The weights
of the probabilities are either 1/3 or 2/3 at the two boundaries, while they
are 1/9, 2/9, or 4/9 at the corners.

Decimating the sampling number of PSFs is an effective approach to minimize
computational cost. However, this technique can introduce side-effects at the
boundary of segments when switching between PSFs. The variation in shape
between neighboring PSFs, caused by the sampling interval, leads to
discontinuities in the deconvolution process. To mitigate this issue, a simple
countermeasure is to randomly select adjacent PSFs at their boundaries, which
helps to smooth out the discontinuities. The weights of the probabilities for
selecting PSFs are illustrated in Figure 7. The presence and severity of
artifacts depend on factors such as the dissimilarity in shape between
neighboring PSFs, statistical characteristics of the observed images, and
other relevant factors. Therefore, the presented technique serves as an
example, and the problem is optimized by the range of pixels to be randomized.

Figure 8: Comparison of the RLsv method without and with PSF randomization at
the boundaries. (a) PSF images corresponding to (b) and (c). (b) RLsv-image
without correcting the PSF boundaries. (b-1) Enlarged images specified by the
colored frames in (b). The arrows correspond to the boundaries of the PSFs.
(c) Same as (b), but using the randomization of the PSFs.

The result of applying the selection rule in Figure 7 is shown in Figure 8. We
compared the RLsv method for the observed data in the eastern region in
Section 3.1 with a PSF sampling of $35\times 35$ pixels and 200 iterations.
Figure 8(a) shows the PSF images corresponding to Figures 8(b) and (c), where
the white lines are used to clarify the border lines. Figures 8(b) and (c) are
the RLsv-deconvolved images with and without using the randomization of PSFs,
respectively. To illustrate the differences, we presented magnified images of
Figures 8(b) and (c) as Figures 8(b-1) and (c-1), respectively. It appears
that the discontinuity at the boundaries of the PSFs, indicated by the green
arrows, is smeared out to some extent.

## References

  * Agarwal et al. (2020) Agarwal, C., Khobahi, S., Bose, A., Soltanalian, M., & Schonfeld, D. 2020, in 2020 IEEE International Conference on Image Processing (ICIP), 3299–3303, doi: 10.1109/ICIP40778.2020.9190825
  * Bamba et al. (2003) Bamba, A., Yamazaki, R., Ueno, M., & Koyama, K. 2003, The Astrophysical Journal, 589, 827
  * Burke et al. (1994) Burke, B., Mountain, R., Daniels, P., Cooper, M., & Dolat, V. 1994, IEEE Transactions on Nuclear Science, 41, 375
  * Dalitz et al. (2015) Dalitz, C., Pohle-Frohlich, R., & Michalk, T. 2015, IEEE transactions on ultrasonics, ferroelectrics, and frequency control, 62, 531
  * Davis et al. (2012) Davis, J. E., Bautz, M. W., Dewey, D., et al. 2012, in Space Telescopes and Instrumentation 2012: Ultraviolet to Gamma Ray, Vol. 8443, SPIE, 375–386
  * Dey et al. (2006) Dey, N., Blanc-Feraud, L., Zimmer, C., et al. 2006, Microscopy research and technique, 69, 260
  * Esch et al. (2004) Esch, D. N., Connors, A., Karovska, M., & van Dyk, D. A. 2004, The Astrophysical Journal, 610, 1213
  * Fabbiano et al. (2020) Fabbiano, G., Paggi, A., Karovska, M., et al. 2020, The Astrophysical Journal, 902, 49
  * Fruscione et al. (2006) Fruscione, A., McDowell, J. C., Allen, G. E., et al. 2006, in Observatory Operations: Strategies, Processes, and Systems, Vol. 6270, SPIE, 586–597
  * Grefenstette et al. (2015) Grefenstette, B. W., Reynolds, S. P., Harrison, F. A., et al. 2015, The Astrophysical Journal, 802, 15
  * Hwang et al. (2004) Hwang, U., Laming, J. M., Badenes, C., et al. 2004, The Astrophysical Journal, 615, L117
  * Jerius et al. (2000) Jerius, D., Donnelly, R. H., Tibbetts, M., et al. 2000, in X-Ray Optics, Instruments, and Missions III, Vol. 4012, SPIE, 17–27
  * Ku (1966) Ku, H. H. 1966, Journal of Research of the National Bureau of Standards, 70, 263
  * Lucy (1974) Lucy, L. B. 1974, The astronomical journal, 79, 745
  * Markevitch et al. (2000) Markevitch, M., Ponman, T., Nulsen, P., et al. 2000, The Astrophysical Journal, 541, 542
  * Naik & Sahu (2013) Naik, R. K., & Sahu, P. 2013, in 2013 International Conference on Microwave and Photonics (ICMAP), IEEE, 1–6
  * Neelamani et al. (2004) Neelamani, R., Choi, H., & Baraniuk, R. 2004, IEEE Transactions on signal processing, 52, 418
  * Paggi et al. (2022) Paggi, A., Fabbiano, G., Nardini, E., et al. 2022, The Astrophysical Journal, 927, 166
  * Patnaude & Fesen (2009) Patnaude, D. J., & Fesen, R. A. 2009, The Astrophysical Journal, 697, 535
  * Richardson (1972) Richardson, W. H. 1972, JoSA, 62, 55
  * Rudin et al. (1992) Rudin, L. I., Osher, S., & Fatemi, E. 1992, Physica D: nonlinear phenomena, 60, 259
  * Sato et al. (2018) Sato, T., Katsuda, S., Morii, M., et al. 2018, The Astrophysical Journal, 853, 46
  * Sekko et al. (1999) Sekko, E., Thomas, G., & Boukrouche, A. 1999, Signal processing, 72, 23
  * Shepp & Vardi (1982) Shepp, L. A., & Vardi, Y. 1982, IEEE transactions on medical imaging, 1, 113
  * Sobolenko et al. (2022) Sobolenko, M., Kompaniiets, O., Berczik, P., et al. 2022, Monthly Notices of the Royal Astronomical Society, 517, 1791
  * Starck et al. (2002) Starck, J.-L., Pantin, E., & Murtagh, F. 2002, Publications of the Astronomical Society of the Pacific, 114, 1051
  * Sugizaki et al. (2009) Sugizaki, M., Kamae, T., & Maeda, Y. 2009, Publications of the Astronomical Society of Japan, 61, S55
  * Tajima et al. (2007) Tajima, H., Finazzi, S., Cohen-Tanugi, J., Chiang, J., & Kamae, T. 2007in , American Institute of Physics, 187–189
  * Thimmappa et al. (2020) Thimmappa, R., Marchenko, V., Balasubramaniam, K., et al. 2020, The Astrophysical Journal, 903, 109
  * Tsuchioka et al. (2022) Tsuchioka, T., Sato, T., Yamada, S., & Uchiyama, Y. 2022, The Astrophysical Journal, 932, 93
  * van Kempen & van Vliet (2000) van Kempen, G. M., & van Vliet, L. J. 2000, JOSA A, 17, 425
  * White (1994) White, R. L. 1994, in Instrumentation in Astronomy VIII, Vol. 2198, SPIE, 1342–1348
  * Wise et al. (1997) Wise, M. W., Huenemoerder, D. P., & Davis, J. E. 1997, in Astronomical Data Analysis Software and Systems VI, Vol. 125, 477
  * Yongpan et al. (2010) Yongpan, W., Huajun, F., Zhihai, X., Qi, L., & Chaoyue, D. 2010, Optics & Laser Technology, 42, 845
  * Yuan et al. (2008) Yuan, L., Sun, J., Quan, L., & Shum, H.-Y. 2008, Acm Transactions on Graphics (TOG), 27, 1


# NLO corrections to $J/\psi+c+\bar{c}$ photoproduction

Qi-Ming Feng1 and Cong-Feng<EMAIL_ADDRESS>1 School of Physics,
University of Chinese Academy of Sciences, Beijing 100049, China
2 Key Laboratory of Vacuum Physics of CAS, Beijing 100049, China

###### Abstract

Based on the factorization framework of nonrelativistic quantum
chromodynamics, we study the associated $J/\psi+c+\bar{c}$ photoproduction
process at next-to-leading order in $\alpha_{s}$ and leading order in the
velocity expansion. The total cross section and differential cross section in
$p_{T}^{2}$, $W$ and $z$ are presented. The results indicate that the next-to-
leading order corrections are substantial, and testable in experiment.

##  I Introduction

The study of heavy quarkonium system presents an exceptional opportunity to
explore the nuances of the phenomena quantum chromodynamics (QCD) involved.
Especially, charmonium provides a peculiar playground for the study of flavor
physics and QCD in the charm sector, of which a huge amount of experimental
data have been accumulated. The moderate charm energy enables the perturbative
QCD(pQCD) calculation reliable to some extent, while poses a challenge to
higher order pQCD corrections.

Nonrelativistic QCD (NRQCD) Bodwin:1994jh  provides a consistent theoretical
framework for the study of quarkonium production and decays. In NRQCD, the
quarkonium production and decay processes are factorized into two sectors: the
perturbative generation and decay of heavy quark pairs, the dominant
quarkonium quark components, and the nonperturbative hadronization or
dehadronization of these heavy quarks. The perturbative contributions are
represented by the matching coefficients to pQCD calculation, namely the
short-distance coefficients (SDCs), while the nonperturbative hadronization is
described by the matrix elements of process-independent effective operators,
known as the long-distance matrix elements (LDMEs).

Nevertheless, there are still some unsolved problems pending in the
application and understanding of NRQCD. The NRQCD factorization formalism
suggests that besides the leading Fock state contribution usually
corresponding to the color singlet mechanism (CSM), higher Fock state
contributions, such as the color octet mechanism (COM), emerge in the
expansion of heavy quark relative velocity $v$ ($v\ll 1$). The proposal of COM
effectively reduces the discrepancies in $J/\psi$ production between next-to-
leading order (NLO) CSM predictions and experimental results across
$e^{+}e^{-}$ collisions at B factories, photoproduction at DESY HERA, and
hadroproduction at Fermilab Tevatron and CERN LHC Chang:2009uj ;
Artoisenet:2009xh ; Campbell:2007ws ; Gong:2008sn ; Lansberg:2010vq ;
Kramer:1994zi ; Kramer:1995nb ; Butenschoen:2009zy ; Butenschoen:2011ks ;
Butenschoen:2012px ; Chao:2012iv ; Ma:2010jj ; Ma:2010yw ; Zhang:2009ym ;
Gong:2012ug . However, the COM introduces considerable uncertainties. In Ref.
Bodwin:2012ft , comparisons between two LDMEs fitted through different
procedures in various collision processes show somehow incompatible results.
Since different processes rely on distinct sets of LDME data, the process-
independence of COM LDMEs is challenged. Some new methods for fitting COM
LDMEs have been proposed later, but discussing them in detail is beyond the
scope of this text and will be skipped here.

Experimental and theoretical inquiries into inclusive heavy quarkonium
production have spanned several decades (see Brambilla:2010cs ;
Lansberg:2019adr ; QuarkoniumWorkingGroup:2004kpm  for reviews). Recent
studies of $J/\psi$ photoproduction in electron-proton ($ep$) collisions
indicate that CS contributions, such as intrinsic charm Flore:2020jau  and
higher-order processes like $J/\psi+c+\bar{c}$ Li:2019nlr , are evident.
Inspired by the fact that production processes akin to NLO QCD corrections
exhibit notable contributions Campbell:2007ws ; Chen:2016hju ; Yang:2022yxb ,
we posit that the NLO contributions of the aforementioned 3-body final states
photoproduction process $\gamma+g\to J/\psi+c+\bar{c}$ at $ep$ colliders
remains relatively significant. Since the $J/\psi+c+\bar{c}$ final state is
experimentally detectable, theoretical analysis of HERA data holds
significance and provides insights for future $ep$ colliders like EIC, EicC,
and LHeC (FCC-eh).

In this work, based on the framework of NRQCD, we systematically compute the
NLO corrections to the photoproduction process $\gamma+g\to J/\psi+c+\bar{c}$
at leading $v$ expansion in $ep$ collisions. The structure of this paper is
organized as follows. In Section II, we detail the formalism and calculation
of the concerned process. In Section III, the results of numerical evaluation
are presented. The last section is reserved for the summary and conclusions.

##  II Formalism and Calculation

Within the framework of NRQCD, the cross section for the photoproduction
process at leading $v$ in $ep$ collision can be formulated as:

$\displaystyle d\sigma(ep\to J/\psi+c+\bar{c})=\int$ $\displaystyle dxd\eta
f_{\gamma/e}(x,Q^{2}_{\max})f_{g/p}(\eta,\mu^{2})$ $\displaystyle\\!\\!\times
d\sigma(\gamma+g\to
c\bar{c}[^{3}\\!S_{1}]+c+\bar{c})\langle\mathcal{O}^{J/\psi}(^{3}\\!S_{1})\rangle\
.$ (1)

Here, $\langle\mathcal{O}^{J/\psi}(^{3}\\!S_{1})\rangle$ is the LDME of a
$c\bar{c}[^{3}\\!S_{1}]$ pair hadronizing into a $J/\psi$ meson.
$f_{g/p}(\eta,\mu^{2})$ is the parton distribution function (PDF) of the
incident gluon, and $\mu$ is the corresponding factorization scale.
$f_{\gamma/e}(x,Q^{2}_{\max})$ is the Weizsacker-Williams approximation (WWA)
function of the photon distribution, defined as:

$\displaystyle
f_{\gamma/e}(x,Q^{2}_{\max})=\frac{\alpha}{2\pi}\left[\frac{1+{(1-x)}^{2}}{x}\ln\frac{Q^{2}_{\max}}{Q^{2}_{\min}(x)}+2m_{e}^{2}x\left(\frac{1}{Q^{2}_{\max}}-\frac{1}{Q^{2}_{\min}(x)}\right)\right]\
,$ (2)

where $m_{e}$ is the electron mass. $Q^{2}_{\min}=m_{e}^{2}x^{2}/(1-x)$ and
$Q^{2}_{\max}$ represents the minimal and maximum virtuality of the incident
photon, respectively.

In the calculation of the concerned process, the spinor helicity formalism
Kleiss:1985yh ; Qiao:2003ue ; Dixon:1996wi ; Dixon:2013uaa ; Arkani-
Hamed:2017jhn  of the scattering amplitude and the conventional amplitude
squaring approach are introduced in evaluating diagrams. Specifically, most of
the diagrams are calculated in helicity amplitudes, except for the Coulomb
divergent part, where the conventional amplitude squaring approach is
employed.

The dipole subtraction method Catani:1996vz ; Catani:2002hc  is adopted to
counteract the infrared (IR) poles. Therefore, the total cross section at NLO
can be expressed as:

$\displaystyle\sigma_{tot}$
$\displaystyle=\int_{3\text{-}\mathrm{body}}(d\sigma^{\mathrm{LO}}+d\sigma^{\mathrm{Virtual}}+d\sigma^{C}+\int_{1\text{-}\mathrm{body}}d\sigma^{A})+\int_{4\text{-}\mathrm{body}}(d\sigma^{\mathrm{Real}}-d\sigma^{A})\
,$ (3)

where $d\sigma^{\mathrm{LO}}$, $d\sigma^{\mathrm{Virtual}}$, and
$d\sigma^{\mathrm{Real}}$ are the LO, virtual, and real contributions to the
cross section. $d\sigma^{C}$ represents the collinear subtraction counterterm
arising from the redefination of the parton distributions. The contribution
$d\sigma^{A}$ represents the dipole counterterm. The correspondence of the
various terms in (3) with helicity amplitudes goes as

$\displaystyle d\sigma^{\mathrm{LO}}\ ,\ d\sigma^{C}\ ,\
\int_{1\text{-}\mathrm{body}}d\sigma^{A}\propto\sum_{\alpha_{1},\alpha_{2},\alpha_{3},\alpha_{4},\alpha_{5}}|\mathcal{A}_{\mathrm{LO}}^{\alpha_{1},\alpha_{2},\alpha_{3},\alpha_{4},\alpha_{5}}|^{2}\
,$ $\displaystyle
d\sigma^{\mathrm{Virtual}}\propto\sum_{\alpha_{1},\alpha_{2},\alpha_{3},\alpha_{4},\alpha_{5}}2\mathrm{Re}\left[{(\mathcal{A}_{\mathrm{LO}}^{\alpha_{1},\alpha_{2},\alpha_{3},\alpha_{4},\alpha_{5}})}^{*}\mathcal{A}_{\mathrm{Virtual}}^{\alpha_{1},\alpha_{2},\alpha_{3},\alpha_{4},\alpha_{5}}\right]\
,$ $\displaystyle
d\sigma^{\mathrm{Real}}\propto\sum_{\alpha_{1},\alpha_{2},\alpha_{3},\alpha_{4},\alpha_{5},\alpha_{6}}|\mathcal{A}_{\mathrm{Real}}^{\alpha_{1},\alpha_{2},\alpha_{3},\alpha_{4},\alpha_{5},\alpha_{6}}|^{2}\
.$ (4)

Here, $\alpha_{1,2,3,4,5,6}$ represent the helicities (polarizations) of on-
shell particles, $\alpha_{1,2,4,5,6}\in{(+,-)}$ denote the helicities of
initial photon, initial gluon, the final charm quark, the charm antiquark, and
the final emitted gluon, $\alpha_{3}\in{(+,0,-)}$ signify the helicities of
$J/\psi$ meson.
$\mathcal{A}_{\mathrm{LO}}^{\alpha_{1},\alpha_{2},\alpha_{3},\alpha_{4},\alpha_{5}}$,
$\mathcal{A}_{\mathrm{Virtual}}^{\alpha_{1},\alpha_{2},\alpha_{3},\alpha_{4},\alpha_{5}}$,
and
$\mathcal{A}_{\mathrm{Real}}^{\alpha_{1},\alpha_{2},\alpha_{3},\alpha_{4},\alpha_{5},\alpha_{6}}$
are helicity amplitudes of LO, virtual corrections, and real corrections,
respectively.

Figure 1: Typical LO Feynman diagrams for $\gamma+g\to J/\psi+c+\bar{c}$. All
other diagrams can be generated by: 1. exchanging initial photon and gluon in
$(a)$, $(b)$ and $(d)$; 2. reversing fermion lines; 3. constraining $J/\psi$
in other quark-antiquark pairs. There are 6 diagrams related to (a), 12
diagrams related to (b), 4 diagrams related to (c), and 8 diagrams related to
(d). Diagrams with no contribution in the end are neglected.

There are 30 non-zero Feynman diagrams at the leading order of the concerned
process, as schematically shown in FIG. 1. Helicity amplitude of the LO
contribution can be expressed as:

$\displaystyle\mathcal{A}_{\mathrm{LO}}^{\alpha_{1},\alpha_{2},\alpha_{3},\alpha_{4},\alpha_{5}}=\sum_{i=1}^{30}\mathcal{A}^{\alpha_{1},\alpha_{2},\alpha_{3},\alpha_{4},\alpha_{5}}_{i}=\sum_{j=1}^{59}\mathcal{C}_{j}(p_{a}\cdot
p_{b},p_{a}\cdot\varepsilon^{\alpha_{b}}_{b},\varepsilon^{\alpha_{a}}_{a}\cdot\varepsilon^{\alpha_{b}}_{b})\mathcal{\hat{A}}_{j}^{\alpha_{1},\alpha_{2},\alpha_{3},\alpha_{4},\alpha_{5}}\
,$ (5)

where $p_{a}$ and $p_{b}$ represent on-shell momenta,
$\varepsilon^{\alpha_{a}}_{a}$ and $\varepsilon^{\alpha_{b}}_{b}$ represent
polarization vectors. The total helicity amplitude
$\mathcal{A}_{\mathrm{LO}}^{\alpha_{1},\alpha_{2},\alpha_{3},\alpha_{4},\alpha_{5}}$
sums over all 30 LO diagrams
$\mathcal{A}_{i}^{\alpha_{1},\alpha_{2},\alpha_{3},\alpha_{4},\alpha_{5}}$. We
rearrange the summation into 59 distinct combinations of spinor products for
simplification, expressed as
$\mathcal{\hat{A}}_{j}^{\alpha_{1},\alpha_{2},\alpha_{3},\alpha_{4},\alpha_{5}}$.
$\mathcal{C}_{j}$ represents corresponding coefficients of
$\mathcal{\hat{A}}_{j}^{\alpha_{1},\alpha_{2},\alpha_{3},\alpha_{4},\alpha_{5}}$,
composed of scalar products of momenta and polarization vectors.

In the calculation of NLO corrections, the ’t Hooft-Veltman scheme (HV)
dimensional regularization (DR) is employed. Some of the NLO diagrams are
shown in FIG. 2.

Figure 2: Typical NLO Feynman diagrams of $\gamma+g\to J/\psi+c+\bar{c}$.
$(a)$ and $(b)$ are typical renormalizaiton diagrams, where the $\otimes$
denotes the counterterm of a propagator or a vertex. $(c)\sim(e)$ are typical
loop diagrams with UV, IR and Coulomb poles, respectively. The UV pole in (c)
can be canceled by introducing renormalization. The IR pole in $(d)$ can be
counteracted with real corrections (initial gluon introduced pole) and other
loop diagrams (the pole introduced form $J/\psi$ and $c$). The coulomb pole in
$(e)$ can be factorization out. And $(f)$ is a typical diagram of the real
correction process $\gamma+g\to J/\psi+c+\bar{c}+g$.

Loop diagrams are evaluated using two distinct methods, the integrand
reduction and IBP reduction. Helicity amplitudes of loop diagrams without
Coulomb divergence are evaluated using integrand reduction via Laurent-
expansion method Mastrolia:2012bu  by the semi-numerical reduction C++ library
Ninja Peraro:2014cba . For the Coulomb divergent diagrams, the Coulomb term
can be expressed as:

$\displaystyle\frac{1}{[{(-\frac{p_{3}}{2}+q)}^{2}-m_{c}^{2}]q^{2}[(\frac{p_{3}}{2}+q)^{2}-m_{c}^{2}]}=\frac{1}{2}\left(\frac{1}{{(q^{2})}^{2}[(-\frac{p_{3}}{2}+q)^{2}-m_{c}^{2}]}+\frac{1}{{(q^{2})}^{2}[{(\frac{p_{3}}{2}+q)}^{2}-m_{c}^{2}]}\right),~{}~{}$
(6)

where only the loop propagators related to the Coulomb poles are depicted.
$p_{3}$ is $J/\psi$ meson momentum, $q$ is the loop momentum, and $m_{c}$ is
the charm quark mass. The introduction of exceptional higher-power loop
propagators $\frac{1}{{(q^{2})}^{2}}$ prevents the evaluation of Coulomb
divergent diagrams using Ninja. Therefore, we employ the integration by parts
(IBP) reduction method to evaluate Coulomb loops by NeatIBP Wu:2023upw .

Ultraviolet (UV) singularities are canceled by renormalization. In our
calculation, renormalizaiton constant of the QCD coupling constant $Z_{g}$ is
defined in the modified minimal subtraction ($\overline{\rm MS}$) scheme,
renormalizaiton constant of the charm quark field $Z_{2}$ and mass $Z_{m}$ and
the gluon field $Z_{3}$ are defined in the on-shell (OS) schem. The
counterterms are given by:

$\displaystyle\delta Z_{2}^{\rm
OS}=-C_{F}\frac{\alpha_{s}}{4\pi}\left[\frac{1}{\epsilon_{\rm{UV}}}+\frac{2}{\epsilon_{\rm{IR}}}-3\gamma_{E}+3\ln\frac{4\pi\mu_{r}^{2}}{m^{2}}+4\right]\
,$ $\displaystyle\delta Z_{3}^{\rm
OS}=\frac{\alpha_{s}}{4\pi}\left[(\beta^{\prime}_{0}-2C_{A})(\frac{1}{\epsilon_{\rm{UV}}}-\frac{1}{\epsilon_{\rm{IR}}})-\frac{4}{3}T_{f}\left(\frac{1}{\epsilon_{\rm{UV}}}-\gamma_{E}+\ln\frac{4\pi\mu_{r}^{2}}{m^{2}}\right)\right]\
,$ $\displaystyle\delta Z_{m}^{\rm
OS}=-3C_{F}\frac{\alpha_{s}}{4\pi}\left[\frac{1}{\epsilon_{\rm
UV}}-\gamma_{E}+\ln\frac{4\pi\mu_{r}^{2}}{m^{2}}+\frac{4}{3}\right]\ ,$
$\displaystyle\delta Z_{g}^{\overline{\rm
MS}}=-\frac{\beta_{0}}{2}\frac{\alpha_{s}}{4\pi}\left[\frac{1}{\epsilon_{\rm
UV}}-\gamma_{E}+\ln(4\pi)\right]\ .$ (7)

Here, $\beta^{\prime}_{0}=(11/3)C_{A}-(4/3)T_{f}n^{\prime}_{f}$ is the one-
loop coefficient of the QCD beta function, $n^{\prime}_{f}=3$ is the number of
light quarks, $\gamma_{E}$ is Euler’s constant, $m$ represents the mass of
charm quark, $C_{A}$ and $T_{f}$ attribute to the color $SU(3)$ group,
$n_{f}=4$ is the number of active quarks, $\mu_{r}$ denotes the
renormalization scale.

In the evaluation of real corrections, as expected, the dipole counterterm
$d\sigma^{A}$ exhibits the same pointwise singular behaviour as
$d\sigma^{\mathrm{Real}}$, and $\int_{1\text{-}\mathrm{body}}d\sigma^{A}$
stands for the analytic integration of $d\sigma^{A}$ over the phase space with
an additional real gluon in dimension $D=4-2\varepsilon$, cancelling out the
remaining analytic $\frac{1}{\varepsilon}$ and $\frac{1}{\varepsilon^{2}}$
divergences in virtual correction.

In the concerned process, the dipole terms associated with quarkonium cancel
out. As a result, only 3 types of dipole terms remain:

  1. 1.

initial gluon emitter with final charm (anti-charm) spectator:
$\mathcal{D}^{gg}_{c},\ \mathcal{D}^{gg}_{\bar{c}}$,

  2. 2.

final charm (anti-charm) emitter with initial gluon spectator:
$\mathcal{D}^{g}_{cg},\ \mathcal{D}^{g}_{\bar{c}g}$,

  3. 3.

final charm (anti-charm) emitter with final anti-charm (charm) spectator:
$\mathcal{D}_{cg,\bar{c}},\ \mathcal{D}_{\bar{c}g,c}$,

Here, dipole contributions $\mathcal{D}^{gg}_{c}$,
$\mathcal{D}^{gg}_{\bar{c}}$, $\mathcal{D}^{g}_{cg}$,
$\mathcal{D}^{g}_{\bar{c}g}$, $\mathcal{D}_{cg,\bar{c}}$, and
$\mathcal{D}_{\bar{c}g,c}$ are defined in Ref. Catani:2002hc . Hence, the
dipole factorization form of $d\sigma^{A}$ writes:

$\displaystyle
d\sigma^{A}=d\Gamma^{(4)}\left(\sum_{\begin{subarray}{c}i,k=c,\bar{c}\\\ i\neq
k\end{subarray}}\mathcal{D}_{ig,k}+\sum_{i=c,\bar{c}}\mathcal{D}_{ig}^{g}+\sum_{k=c,\bar{c}}\mathcal{D}_{k}^{gg}\right)$
(8)

with $d\Gamma^{(4)}$ being the 4-body phase space including all factors of QCD
independent. The integrated dipoles are obtained from $(5.23)$, $(5.56)$, and
$(5.88)$ of Ref. Catani:2002hc .

Of the cancellation of divergences, we extract the IR divergences in loop
diagrams by means of the method developed in Ref. Dittmaier:2003bc , which
ensures all IR divergences in one-loop diagrams with more than three loop
propagators being expressed as sum of triangles (diagrams with three loop
propagators). For example, for an $N$-point loop ($N\geq 3$)

$\displaystyle\int_{\rm{div}}d^{4-2\epsilon}q\prod_{i=0}^{N-1}\frac{1}{D_{i}}=\int_{\rm{div}}d^{4-2\epsilon}q\sum_{i=0}^{N-1}\sum_{\begin{subarray}{c}j=0\\\
k\neq i,i+1\end{subarray}}^{N-1}\frac{A_{ij}}{D_{i}D_{i+1}D_{j}}\ ,$ (9)

where $D_{i}$ represents loop propagators. $A_{ij}$ represents the
corresponding coefficients. In the calculation, divergences in virtual
correction and dipole contribution are analytically canceled out.

##  III Results

In our numerical calculation, the charm mass takes half of $J/\psi$ mass, i.e.
$m_{c}=1.5\ \rm{GeV}$. The renormalization scale $\mu_{r}$ and the
factorizaiton scale $\mu_{f}$ are set as
$\mu_{r}=\mu_{f}=m_{T}\equiv\sqrt{p_{T}^{2}+4m_{c}^{2}}$, where $m_{T}$ is the
$J/\psi$ transverse mass. Theoretical uncertainties are estimated by varying
the charm quark mass in $m_{c}=1.5\pm 0.1\ \rm{GeV}$ and the scales in the
interval $\frac{1}{2}m_{T}\leq\mu_{r},\mu_{f}\leq 2m_{T}$. The running
coupling constant is determined by the one-loop (two-loop) formula at LO
(NLO), and the PDF set CTEQ6L1 (CTEQ6M) Pumplin:2002vw  is used at LO (NLO).
The LDME follows
$\langle\mathcal{O}^{J/\psi}(^{3}\\!S_{1})\rangle=2(2J+1)N_{c}|R(0)|^{2}/{4\pi}$
with $J=1$ for the $J/\psi$ meson, and $N_{c}=3$ is the number of color
charges. The radial wave function $|R(0)|^{2}=1.01\ \rm{GeV}^{3}$ is extract
form $\Gamma(J/\psi\to e^{+}e^{-})=5.55\ \rm{keV}$ ParticleDataGroup:2020ssz ,
thus, we have $\langle\mathcal{O}^{J/\psi}(^{3}\\!S_{1})\rangle=1.45\
\rm{GeV}^{3}$.

The collision energy is set accroding to the HERA collider: $27.5\ \rm{GeV}$
for eletrons (positrons) and $920\ \rm{GeV}$ for protons. The photon
virtuality is constrainted to $Q^{2}_{\max}\leq 2.5\ \rm{GeV}^{2}$. To exclude
resolved photoproduction and diffractive production of $J/\psi$, experimental
cuts based on the H1 collaboration measurement H1:2010udv  are applied:
$p_{T}>1\ \rm{GeV}$, $60\ \rm{GeV}<W<240\ \rm{GeV}$, and $0.3<z<0.9$. Here,
$W=\sqrt{{(p_{\gamma}+p_{p})}^{2}}$ is the mass of the hadronic final state,
$z=(p_{3}\cdot p_{p})/(p_{\gamma}\cdot p_{p})$ is the elasticity of the
$J/\psi$ meson production process, and $p_{p}$, $p_{\gamma}$ are the momenta
of the incident proton and photon, respectively. Additionally, the feed-down
contribution from the $\psi^{\prime}$ is taken into account, which yields an
enhancing factor about 0.278.

Table 1: Scale and mass dependence of the total cross section at LO (expressed in $\rm{nb}$) in various PDF sets without feed-down contribution. Here, $\mu=\mu_{r}=\mu_{f}$.  PDF sets |  $m_{c}\backslash\mu$ |  $\frac{1}{2}m_{T}$ |  $m_{T}$ |  $2m_{T}$
---|---|---|---|---
CTEQ6L1 |  $1.4\ \rm{GeV}$ |  $0.185$ |  $0.111$ |  $0.070$
|  $1.5\ \rm{GeV}$ |  $0.123$ |  $0.074$ |  $0.046$
|  $1.6\ \rm{GeV}$ |  $0.084$ |  $0.049$ |  $0.031$
CT14LO |  $1.4\ \rm{GeV}$ |  $0.128$ |  $0.073$ |  $0.046$
|  $1.5\ \rm{GeV}$ |  $0.084$ |  $0.048$ |  $0.030$
|  $1.6\ \rm{GeV}$ |  $0.056$ |  $0.032$ |  $0.020$
CTEQ6M |  $1.4\ \rm{GeV}$ |  $0.109$ |  $0.066$ |  $0.043$
|  $1.5\ \rm{GeV}$ |  $0.073$ |  $0.044$ |  $0.029$
|  $1.6\ \rm{GeV}$ |  $0.051$ |  $0.030$ |  $0.019$
CT14NLO |  $1.4\ \rm{GeV}$ |  $0.095$ |  $0.063$ |  $0.042$
|  $1.5\ \rm{GeV}$ |  $0.066$ |  $0.042$ |  $0.028$
|  $1.6\ \rm{GeV}$ |  $0.045$ |  $0.029$ |  $0.019$
CT18NLO |  $1.4\ \rm{GeV}$ |  $0.094$ |  $0.063$ |  $0.042$
|  $1.5\ \rm{GeV}$ |  $0.065$ |  $0.042$ |  $0.028$
|  $1.6\ \rm{GeV}$ |  $0.045$ |  $0.029$ |  $0.019$

As a result, the total cross section of the concerned process at NLO (LO) is

$\displaystyle\sigma_{tot}=0.118^{+0.168}_{-0.065}\ (0.074^{+0.111}_{-0.043})\
\rm{nb}.$ (10)

Our LO result agrees with what in Ref. Li:2009zzu  after taking the same
inputs. NLO corrections yield a $K$ factor about $1.60$, which is a prominent
enhancement of the cross section. It is evident that the error on the total
cross section is large. In Table 1, we present the errors of the LO total
cross section on the charm quark mass and scales across various PDF sets. The
cross section with each PDF set shows strong sensitivity to both charm quark
mass and scales, particularly to scales. Results using the NLO PDF sets
CTEQ6M, CT14NLO Dulat:2015mca , and CT18NLO Hou:2019qau  exhibit stronger
dependence on scales compared to those obtained using LO PDF sets CTEQ6L1 and
CT14LO Dulat:2015mca .

In estimating the number of the concerned process events at HERA, we consider
that the H1 collaboration reconstructed the $J/\psi$ meson candidates through
the decay channel $J/\psi\to\mu^{+}\mu^{-}$ in the photoproduction process,
and the photoproduction sample corresponds to an integrated luminosity of
$\mathcal{L}=165\ \rm{pb}^{-1}$ H1:2010udv . Accroding to our numerical
results, with a branching fraction of
$\Gamma(J/\psi\to\mu^{+}\mu^{-})/\Gamma_{tot}\simeq 6\%$, the number of
reconstructed $J/\psi+c+\bar{c}$ events in the photoproduction process at NLO
(LO) is about $521\sim 2833\ (304\sim 1829)$ at HERA. Here, we omit the
unclear tagging efficiency of $c$-jets.

Figure 3: The differential cross section in (a) $p_{T}^{2}$, (b) $W$, and (c)
$z$ distributions of the photoproduction process $\gamma+g\to
J/\psi+c+\bar{c}$ at LO and NLO in the CSM. The shaded bands indicate the
theoretical uncertainties with upper bound for
$\mu_{r},\mu_{f}=\frac{1}{2}m_{T}$ and $m_{c}=1.4\ \rm{GeV}$ and lower bound
for $\mu_{r},\mu_{f}=2m_{T}$ and $m_{c}=1.6\ \rm{GeV}$.

The differential cross section distributions in $p_{T}^{2}$, $W$, and $z$, are
presented in FIG. 3 $(a)\sim(c)$, respectively. The $p_{T}^{2}$ distribution
of $J/\psi$, as shown in FIG. 3$(a)$, is presented in the range $1\
\rm{GeV}^{2}<p_{T}^{2}<100\ \rm{GeV}^{2}$. Compared to the $p_{T}^{2}$
distribution of the inclusive process $\gamma+g\to J/\psi+X$ at NLO CSM in
Ref. Butenschoen:2009zy , the $J/\psi+c+\bar{c}$ process have a much milder
drop as $p_{T}^{2}$ increases. Although the differential cross section of the
concerned process is much lower than that of the $J/\psi+X$ production process
at low-$p_{T}$, in the region of $60\ \rm{GeV}<p_{T}^{2}<100\ \rm{GeV}^{2}$,
the ratio $\sigma(\gamma+g\to J/\psi+c+\bar{c})/\sigma(\gamma+g\to J/\psi+X)$
tends to about $1/2$. It means that the concerned process is a not negligible
contribution in the evaluation of $J/\psi$ photoproduction processes in $ep$
colliders, particularly in the large-$p_{T}$ region.

The $W$ and $z$ distributions in FIG. 3$(b)\sim(c)$ significantly undershoot
the inclusive $\gamma+g\to J/\psi+X$ process in the CSM in Ref.
Butenschoen:2009zy . The $W$ distribution of the concerned process, the CS
$J/\psi+X$ production process, and the H1 measurement show a similar trend.
However, the $z$ distribution shows a completely different trend compared to
the CS $J/\psi+X$ production process and H1 data. In contrast to the moderate
$z$ distribution trends in Ref.Butenschoen:2009zy  and Ref.H1:2010udv , the
concerned process exhibits a rapid decrease in distribution with increasing
$z$, particularly at large $z$.

The discrepancy on $z$ distribution may originate from the distinct dynamics
of these processes. With the elasticity $z$ increases, the $J/\psi$ meson
tends to parallel with the incident photon, that makes the momentum of the sum
of the final charm quark and final anti-charm quark reduced. Since the mass of
the sum momentum greater than the sum of charm quark mass and anti-charm quark
mass: $m_{c+\bar{c}}=\sqrt{{(p_{c}+p_{\bar{c}})}^{2}}>2m_{c}$, the dynamics of
the two (anti-) charm quarks are constrainted, apparently, the cross section
is also limited. For the $\gamma+g\to J/\psi+X$ process, $X=g,q$ is a light
particle, there is no such constraint in $z$, so the $z$-distribution is flat.
This can be confirmed by comparing upper, middle, and lower bounds. As $m_{c}$
drops from $1.6\ \rm{GeV}$ to $1.4\ \rm{GeV}$, the distribution tends to be
more flat.

##  IV Summary

In this work, we investigate the photoproduction of $J/\psi+c+\bar{c}$ in $ep$
collider HERA at NLO QCD and LO in $v$ within the framework of NRQCD.
Numerical result shows that the NLO QCD corrections are prominent in the
concerned process. We also present the differential cross section in $p_{T}$,
$W$ and $z$. The concerned process shows a significantly contribution in
$J/\psi$ photoproduction processes at large-$p_{T}$. Comparing with the
$p_{T}$ distribution of the inclusive $\gamma+g\to J/\psi+X$ process in the
CSM, the concerned process have a much milder drop as $p_{T}$ increases. The
trend in the $W$ distribution typically agrees with that of the CS $J/\psi+X$
photoproduction process and H1 data, by a factor of approximately $1/10$. The
$z$ distribution shows a severe drop, which differs from that of the inclusive
process and H1 data. This difference may originate from the distinct dynamics
of these processes, and expects a further data analysis. Future $ep$
colliders, such as EIC, EicC, and LHeC (FCC-eh), are expected to have rich
production of the concerned process due to their high luminosity. Furthermore,
some of them are capable of detecting particle polarizations, theoretical
analysis on the polarizations of the concerned process awaits future
investigation.

Acknowledgments

This work was supported in part by the National Key Research and Development
Program of China under Contracts Nos.2020YFA0406400, and the National Natural
Science Foundation of China (NSFC) under the Grant 12235008.

## References

  * (1) G. T. Bodwin, E. Braaten and G. P. Lepage, Phys. Rev. D 51 (1995), 1125-1171 [erratum: Phys. Rev. D 55 (1997), 5853] doi:10.1103/PhysRevD.55.5853 [arXiv:hep-ph/9407339 [hep-ph]].
  * (2) C. H. Chang, R. Li and J. X. Wang, Phys. Rev. D 80, 034020 (2009) doi:10.1103/PhysRevD.80.034020 [arXiv:0901.4749 [hep-ph]].
  * (3) P. Artoisenet, J. M. Campbell, F. Maltoni and F. Tramontano, Phys. Rev. Lett. 102, 142001 (2009) doi:10.1103/PhysRevLett.102.142001 [arXiv:0901.4352 [hep-ph]].
  * (4) J. M. Campbell, F. Maltoni and F. Tramontano, Phys. Rev. Lett. 98, 252002 (2007) doi:10.1103/PhysRevLett.98.252002 [arXiv:hep-ph/0703113 [hep-ph]].
  * (5) B. Gong and J. X. Wang, Phys. Rev. Lett. 100, 232001 (2008) doi:10.1103/PhysRevLett.100.232001 [arXiv:0802.3727 [hep-ph]].
  * (6) J. P. Lansberg, Phys. Lett. B 695, 149-156 (2011) doi:10.1016/j.physletb.2010.10.054 [arXiv:1003.4319 [hep-ph]].
  * (7) M. Kramer, J. Zunft, J. Steegborn and P. M. Zerwas, Phys. Lett. B 348, 657-664 (1995) doi:10.1016/0370-2693(95)00155-E [arXiv:hep-ph/9411372 [hep-ph]].
  * (8) M. Krämer, Nucl. Phys. B 459, 3-50 (1996) doi:10.1016/0550-3213(95)00568-4 [arXiv:hep-ph/9508409 [hep-ph]].
  * (9) M. Butenschoen and B. A. Kniehl, Phys. Rev. Lett. 104, 072001 (2010) doi:10.1103/PhysRevLett.104.072001 [arXiv:0909.2798 [hep-ph]].
  * (10) M. Butenschoen and B. A. Kniehl, Phys. Rev. Lett. 108, 172002 (2012) doi:10.1103/PhysRevLett.108.172002 [arXiv:1201.1872 [hep-ph]].
  * (11) M. Butenschoen and B. A. Kniehl, Phys. Rev. Lett. 107, 232001 (2011) doi:10.1103/PhysRevLett.107.232001 [arXiv:1109.1476 [hep-ph]].
  * (12) K. T. Chao, Y. Q. Ma, H. S. Shao, K. Wang and Y. J. Zhang, Phys. Rev. Lett. 108, 242004 (2012) doi:10.1103/PhysRevLett.108.242004 [arXiv:1201.2675 [hep-ph]].
  * (13) Y. Q. Ma, K. Wang and K. T. Chao, Phys. Rev. D 84, 114001 (2011) doi:10.1103/PhysRevD.84.114001 [arXiv:1012.1030 [hep-ph]].
  * (14) Y. Q. Ma, K. Wang and K. T. Chao, Phys. Rev. Lett. 106, 042002 (2011) doi:10.1103/PhysRevLett.106.042002 [arXiv:1009.3655 [hep-ph]].
  * (15) Y. J. Zhang, Y. Q. Ma, K. Wang and K. T. Chao, Phys. Rev. D 81, 034015 (2010) doi:10.1103/PhysRevD.81.034015 [arXiv:0911.2166 [hep-ph]].
  * (16) B. Gong, L. P. Wan, J. X. Wang and H. F. Zhang, Phys. Rev. Lett. 110, no.4, 042002 (2013) doi:10.1103/PhysRevLett.110.042002 [arXiv:1205.6682 [hep-ph]].
  * (17) G. T. Bodwin, [arXiv:1208.5506 [hep-ph]].
  * (18) N. Brambilla, S. Eidelman, B. K. Heltsley, R. Vogt, G. T. Bodwin, E. Eichten, A. D. Frawley, A. B. Meyer, R. E. Mitchell and V. Papadimitriou, et al. Eur. Phys. J. C 71, 1534 (2011) doi:10.1140/epjc/s10052-010-1534-9 [arXiv:1010.5827 [hep-ph]].
  * (19) J. P. Lansberg, Phys. Rept. 889, 1-106 (2020) doi:10.1016/j.physrep.2020.08.007 [arXiv:1903.09185 [hep-ph]].
  * (20) N. Brambilla et al. [Quarkonium Working Group], doi:10.5170/CERN-2005-005 [arXiv:hep-ph/0412158 [hep-ph]].
  * (21) C. Flore, J. P. Lansberg, H. S. Shao and Y. Yedelkina, Phys. Lett. B 811, 135926 (2020) doi:10.1016/j.physletb.2020.135926 [arXiv:2009.08264 [hep-ph]].
  * (22) R. Li, [arXiv:1912.12822 [hep-ph]].
  * (23) Z. Q. Chen, L. B. Chen and C. F. Qiao, Phys. Rev. D 95, no.3, 036001 (2017) doi:10.1103/PhysRevD.95.036001 [arXiv:1608.06231 [hep-ph]].
  * (24) H. Yang, Z. Q. Chen and C. F. Qiao, Phys. Rev. D 105, no.9, 094014 (2022) doi:10.1103/PhysRevD.105.094014 [arXiv:2203.14204 [hep-ph]].
  * (25) R. Kleiss and W. J. Stirling, Nucl. Phys. B 262, 235-262 (1985) doi:10.1016/0550-3213(85)90285-8
  * (26) C. F. Qiao, Phys. Rev. D 67, 097503 (2003) doi:10.1103/PhysRevD.67.097503 [arXiv:hep-ph/0302128 [hep-ph]].
  * (27) L. J. Dixon, [arXiv:hep-ph/9601359 [hep-ph]].
  * (28) L. J. Dixon, doi:10.5170/CERN-2014-008.31 [arXiv:1310.5353 [hep-ph]].
  * (29) N. Arkani-Hamed, T. C. Huang and Y. t. Huang, JHEP 11, 070 (2021) doi:10.1007/JHEP11(2021)070 [arXiv:1709.04891 [hep-th]].
  * (30) S. Catani and M. H. Seymour, Nucl. Phys. B 485, 291-419 (1997) [erratum: Nucl. Phys. B 510, 503-504 (1998)] doi:10.1016/S0550-3213(96)00589-5 [arXiv:hep-ph/9605323 [hep-ph]].
  * (31) S. Catani, S. Dittmaier, M. H. Seymour and Z. Trocsanyi, Nucl. Phys. B 627, 189-265 (2002) doi:10.1016/S0550-3213(02)00098-6 [arXiv:hep-ph/0201036 [hep-ph]].
  * (32) P. Mastrolia, E. Mirabella and T. Peraro, JHEP 06, 095 (2012) [erratum: JHEP 11, 128 (2012)] doi:10.1007/JHEP11(2012)128 [arXiv:1203.0291 [hep-ph]].
  * (33) T. Peraro, Comput. Phys. Commun. 185, 2771-2797 (2014) doi:10.1016/j.cpc.2014.06.017 [arXiv:1403.1229 [hep-ph]].
  * (34) Z. Wu, J. Boehm, R. Ma, H. Xu and Y. Zhang, Comput. Phys. Commun. 295, 108999 (2024) doi:10.1016/j.cpc.2023.108999 [arXiv:2305.08783 [hep-ph]].
  * (35) S. Dittmaier, Nucl. Phys. B 675, 447-466 (2003) doi:10.1016/j.nuclphysb.2003.10.003 [arXiv:hep-ph/0308246 [hep-ph]].
  * (36) J. Pumplin, D. R. Stump, J. Huston, H. L. Lai, P. M. Nadolsky and W. K. Tung, JHEP 07, 012 (2002) doi:10.1088/1126-6708/2002/07/012 [arXiv:hep-ph/0201195 [hep-ph]].
  * (37) P. A. Zyla et al. [Particle Data Group], PTEP 2020, no.8, 083C01 (2020) doi:10.1093/ptep/ptaa104
  * (38) F. D. Aaron et al. [H1], Eur. Phys. J. C 68, 401-420 (2010) doi:10.1140/epjc/s10052-010-1376-5 [arXiv:1002.0234 [hep-ex]].
  * (39) R. Li and K. T. Chao, Phys. Rev. D 79, 114020 (2009) doi:10.1103/PhysRevD.79.114020 [arXiv:0904.1643 [hep-ph]].
  * (40) S. Dulat, T. J. Hou, J. Gao, M. Guzzi, J. Huston, P. Nadolsky, J. Pumplin, C. Schmidt, D. Stump and C. P. Yuan, Phys. Rev. D 93, no.3, 033006 (2016) doi:10.1103/PhysRevD.93.033006 [arXiv:1506.07443 [hep-ph]].
  * (41) T. J. Hou, K. Xie, J. Gao, S. Dulat, M. Guzzi, T. J. Hobbs, J. Huston, P. Nadolsky, J. Pumplin and C. Schmidt, et al. [arXiv:1908.11394 [hep-ph]].