# Accurate Cutting-point Estimation for Robotic Lychee Harvesting through
Geometry-aware Learning

Gengming Zhang, Hao Cao, Kewei Hu, Yaoqiang Pan, Yuqin Deng,
Hongjun Wang#, and Hanwen Kang#
College of Engineering, South China Agriculture University, Guangzhou, China
# Correspond Authors

###### Abstract

Accurately identifying lychee-picking points in unstructured orchard
environments and obtaining their coordinate locations is critical to the
success of lychee-picking robots. However, traditional two-dimensional (2D)
image-based object detection methods often struggle due to the complex
geometric structures of branches, leaves and fruits, leading to incorrect
determination of lychee picking points. In this study, we propose a
Fcaf3d-lychee network model specifically designed for the accurate
localisation of lychee picking points. Point cloud data of lychee picking
points in natural environments are acquired using Microsoft’s Azure Kinect DK
time-of-flight (TOF) camera through multi-view stitching. We augment the Fully
Convolutional Anchor-Free 3D Object Detection (Fcaf3d) model with a squeeze-
and-excitation (SE) module, which exploits human visual attention mechanisms
for improved feature extraction of lychee picking points. The trained network
model is evaluated on a test set of lychee-picking locations and achieves an
impressive $F_{1}$ score of 88.57%, significantly outperforming existing
models. Subsequent three-dimensional (3D) position detection of picking points
in real lychee orchard environments yields high accuracy, even under varying
degrees of occlusion. Localisation errors of lychee picking points are within
±1.5 cm in all directions, demonstrating the robustness and generality of the
model.

##  I Introduction

Precision agriculture or precision farming is a modern concept of whole-farm
management that uses a variety of technologies ranging from remote sensing and
proximal data collection to automation and robotics. In this study we are
looking at the harvesting of fresh lychees. This work has proved to be very
challenging, not only in terms of the mechanical design of the robot, but also
in terms of the vision system, navigation, control methods and manipulation
system [1]. In this case, our main focus is the detection of picking points on
lychee stems in point cloud information to improve the autonomous harvesting
of traditional picking robots.

Single fruits such as apples and oranges can be picked directly, while cluster
fruits such as lychee and grapes need to be picked for the whole cluster.
Lychee can only be picked by identifying the main fruit-bearing branch (MFBB)
picking point and then cutting the MFBB to prevent fruit damage [2]. If the
workspace is clear and fruits are not blocked by obstacles, it is not
challenging to detach these fruits off trees [3]. However, the field is
unstructured, and lychee fruit bunches can be of different sizes in shape and
appear at different heights and positions, in addition, the lychee picking
points in unstructured scenarios are prone to dense occlusion (or the target
position is not obvious enough), and small targets are difficult to detect,
thus affecting the detection accuracy. Thus, it has become the focus of
research to find accurate and directed locating of lychee fruit bunch picking
points, and two key technologies are needed to achieve this demand: (1) three-
dimensional (3D) point cloud data that resist occlusion, and (2) a 3D point
cloud detection CNN model with the ability of precise location under the
premise of the former.

In recent years, computer vision-based technologies have been widely used to
locate picking points for fruit [4, 5, 6]. The characteristics of the target
object itself (such as color, shape, texture, etc.) are used for lychee
detection and identification through image processing (e.g. image filtering,
image segmentation, morphological processing) and machine learning algorithms
and such operations [7, 8, 9, 10, 11]. With the profound development of deep
learning, especially convolutional neural networks (CNN), which further
enhances the generalisation ability as well as the robustness of machine
vision perception systems, more and more researchers have turned to 2D image
detection methods based on deep learning detection. Object detection
algorithms, including You Only Look Once (YOLO) and deep learning networks
improved on YOLO [12, 13, 14, 15, 16], and a two-stage detection represented
by the RCNN (region-based convolutional neural networks) series [17]. The
above study achieved good results for the detection of picking points, but it
is still difficult to challenge the situation due to light variations,
occlusion by leaves and branches, etc.

The advantages of 3D point cloud data itself (e.g. accuracy in representing
geometric shapes, easy access to depth information, 3D reconstruction helps to
improve robustness and accuracy of detection) as well as the intensive
development of 3D convolutional neural networks have motivated researchers to
explore point cloud-based methods for fruit harvesting and picking point
detection. Based on the idea that point cloud reconstruction is not limited by
perspective and helps to improve the robustness and accuracy of detection,
pure manipulation of point cloud data has also become a popular method for
fruit picking point detection [18, 19, 20, 21, 22, 23, 24, 25]. Although the
above research methods have achieved satisfactory results, there are still
significant research gaps in lychee picking point detection.

In this work, a robotic arm-assisted data processing based on TOF camera point
cloud fusion and Fcaf3d-lychee network model for lychee picking points was
proposed, which can be used to accurately locate lychee picking points in the
natural orchard environment, and was tested in the field. The proposed method
fused multi-view point cloud data from a TOF camera, and then visually
perceived the fused sensor data to find the exact location of the lychee
picking point. The results show that the results of this study have higher
localisation accuracy than the typical 3D point cloud target detection model.
In summary, the following contributions are presented in this paper:

  * •

A model Fcaf3d-lychee is proposed for lychee picking point detection, which
greatly improves the accuracy of lychee picking point localization.

  * •

TOF camera for multi-view point cloud stitching, effectively solving the
problem of single-view occlusion and too little sensory information.

  * •

Demonstration of the lychee picking robot equipped with the Fcaf3d-lychee
model on recognition and localization of picking points in the natural orchard
environments.

The rest of this paper is organised as follows. Section II surveys related
work. The system overview and methodologies of our approach are presented in
Section III. The experiment results and discussion are presented in Section
IV, followed by the conclusion in Section V.

##  II Related Works

###  II-A Review on Image-Based 2D Target Detection in lychees

Currently, the application field of robotics is extensive, covering many basic
and related technology areas [26]. Among them, in the intelligent development
of harvesting robots, vision algorithms have been a key factor affecting the
overall performance of the system. The development of vision algorithms for
fruit harvesting robots has two main tasks: target localization as well as
picking point extraction, and the challenges include target color distortion
caused by variable natural light environments, interference of interlaced
plant organs on recognition algorithms, and intraclass deformations of fruits
of the same species [27]. The use of machine vision and its related algorithms
has targeted solutions to improve efficiency, accuracy, intelligence, and
remote interaction during robot operations [28]. To categorize them, the early
major ones are single-feature analysis methods, multi-feature fusion analysis
methods, and pattern recognition methods [29, 30]. [31] proposed a recognition
method based on an improved linear discriminant analysis (LDA) classifier to
address the problem of the low success rate of lychee-picking robots in
recognizing green lychee (not fully ripe lychee) due to background
interference. The core is an LDA method for extracting image convolutional
features, and the idea of ”maximum margin” in the SVM algorithm is introduced
to determine the threshold value of LDA, which is finally integrated into a
multiple LDA classification system by the Adaboost method. Through
experiments, the accuracy of this classifier for unripe lychee is 80.4%, and
the algorithm can also be applied to the determination of fruit ripeness. The
above research is based on image processing, which has high requirements for
light and still needs to go through a certain calculation process when dealing
with the 3D position of the picking point, resulting in poor picking
efficiency and poor robustness.

###  II-B Review on Deep Learning-Based lychee 2D Target Detection Methods

In the existing research, stereo vision measurement techniques and traditional
digital image processing techniques have been adopted by many researchers,
while the booming development of deep learning in the field of image
processing provides a better solution for fruit recognition, thus enabling
fruit-picking robots to adapt to complex orchard environments [32]. To improve
the ability of the picking robot to operate throughout the day, an algorithm
for the detection of lychee fruits and fruit stalks for nighttime environments
was proposed [33]. The performance of the algorithm under different artificial
light intensities was tested separately, and the YOLOv3 convolutional neural
network and U-Net network were utilized to detect fruits and fruit stalks, and
the model finally performed best under normal brightness, achieving an average
detection accuracy of 99.57% as well as segmentation effect of 84.33% MIoU.
However, the method first uses YOLOv3 to obtain the location of the stringed
fruit and then utilizes the ROI region for fruit stalk recognition, which is
not efficient to execute. [34] improved the network structure of YOLOv5s for
the background occlusion and fruit occlusion scenarios that are easily
encountered when picking apples so that it can autonomously determine the
objects that can be grasped and those that cannot be grasped during the
picking operation. The method, however, is mostly applicable to single fruits,
or for estimating yields, and do not apply to lychee harvesting scenarios. The
above studies have made great progress in the acquisition of fruit-picking
points by target detection through 2D images, but most of the studies have not
investigated the direct detection of mother branch picking points.

###  II-C Review on Deep Learning-Based 3D Target Detection Methods for Fruits

Although target detection via 2D images has made great progress in picking
points for fruit-picking robots, most of the studies have not investigated the
direct detection of mother branch picking points. While 2D images may lack
detailed spatial distribution information, 3D point clouds obtained by depth
camera methods provide a more comprehensive representation of object
distribution, especially in complex scenes with occlusions and overlaps. The
potential of 3D point clouds has been increasingly recognized, resulting in
new methods for object detection. [35] proposed a strategy for accurate
picking point localization of artificially cultivated grapes by combining near
and far and utilizing a combination of depth point cloud data. Firstly, the
grape clusters were roughly detected for localization in the far-view image,
and then in the near-view image, the data features of the depth point cloud
and the structural features of the horizontal lattice grapes were fully
utilized, and 95 out of 100 samples were successfully localized, with an
accuracy of 95%. The localization accuracy was as high as 95%. In addition,
[36] proposed a pomegranate tree organ classification and fruit counting
method based on multi-feature fusion and support vector machine (SVM) to
address the limitation of occlusion in the traditional 2D image-based
recognition method, firstly, a 3D point cloud is obtained, then the point
cloud is preprocessed, and its color and shape features are extracted for
training, and a classification model is obtained, and experiments show that
this method can detect most of the fruits on the tree. Existing point cloud
target detection methods, however, have not been studied for picking points.

Our work is based on the Fcaf3d model, which is improved to get the 3D
detection model Fcafd-lychee for lychee picking points, which is able to
utilize the depth camera point cloud stream to directly detect the lychee
picking points and get their 3D coordinates, thus accelerating and efficiently
carrying out the picking point detection.

##  III Materials and Methods

###  III-A System overview

The proposed lychee picking point detection process is shown in Fig. 1. The
process can be divided into two steps, acquiring point cloud and performing
multi-view stitching, and Fcaf3d-lychee modeling for precise location. First,
the picking robot, which has carried the orchard map information, autonomously
walks in the orchard and will obtain the initial position of the target from
the point cloud data stream. After that, the picking robot will perform the
acquisition of the set three-angle point cloud at a close distance to the
target and perform stitching and downsampling. The pre-processed point cloud
will be fed into the Fcaf3d-lychee model for accurate identification and
localization. Finally, the end-effector will be guided to perform the lychee
picking action.

Figure 1: Lychee picking point robot system and methods.

###  III-B Hand-Eye vision model & Closed-Loop calibration method

The camera sensor was operated on the aubo robotic arm with six degrees of
freedom and collected data, which ultimately transformed the 3D coordinates of
the lychee picking point accurately to the coordinate system of the robotic
arm base. Prior to data acquisition, camera calibration and hand-eye
calibration were completed respectively. The former is realized in this study
by the classical Zhang’s calibration method proposed by [37]. In our early
work some research has been done on the calibration between the camera and the
robotic arm, and [38] proposed a robust and accurate closed-loop hand-eye
calibration method to determine the coordinate correlation between the end
flange of the robotic arm and the camera sensor for a robot with an eye in the
hand. The Hand-eye matrix as eq. (1).

$\displaystyle\widehat{{{}_{C}^{F}}T}=\frac{1}{N_{c}}\sum_{i=1}^{N_{C}}{({{}_{B}^{C}}T{{}^{(i)}}{{}_{R}^{B}}T{{}_{F}^{R}}T{{}^{(i)}})^{-1}}$
(1)

In which, the ${}_{R}^{B}T$ is the fixed transformation matrix {R} with
respect to {B} since the relative fixed between the base of the robot arm and
the checkerboard during the calibration process, which can be calibrated by
controlled proximity of the origin of {G} to the origin as well as the X and Y
axes of {B}. the ${}_{B}^{C}T^{(i)}$ is the transformation matrix {B} with
respect to {C}, corresponding to the i th position of the robot that can be
obtained by the least squares method [39]. ${}_{F}^{R}T^{(i)}$ is the
transformation matrix F with respect to R, corresponding to i th position of
the robot that can be easily obtained by the forward kinematics. The
$N{{}_{c}}$ is the number of different poses. In this study, $N{{}_{c}}$=16.

As the robotic arm moved at certain positions [40], the local features at each
location of the lychee picking point were eventually transformed into the same
base coordinate system to make up the global feature (though it should still
be called “local feature” because the arm with the camera sensor was limited
to a localized area when scanning).

###  III-C Point cloud acquisition

####  III-C1 Filtering

Many factors may influence the accuracy of raw point clouds obtained multi-
view vision system, such as changes in natural illumination, vibration in
dynamic environments, hand-eye calibration error, and the camera’s inherent
hardware error. These factors generate messy noise and uneven structures among
the acquired point cloud. As discussed in this section, we combined a
statistical filter and color filter to clear the outliers scattered outside
the main body and with wide color gap from the picking point to provide a good
initial state for subsequent stitching.

As proposed by [41] for filtering for color in 3D reconstruction of fruit
trees, we propose to perform simple color filters in point clouds with color
features to initially reduce a large number of green leaves and provide good
initial conditions for subsequent point cloud detection.

$\displaystyle\left\\{\begin{aligned} R_{s}&>&\sigma_{1}\\\
G_{s}&\leq&\sigma_{2}\end{aligned}\right.$ (2)

where, $R_{s}$, $G_{s}$ are the red and green channels in the color features
of the point cloud respectively, when $R_{s}\leq\sigma_{1}$,
$G_{s}$>$\sigma_{2}$, the point cloud will be removed out of the whole point
cloud.

The statistical filter can effectively remove any anomalous points that are
scattered outside of the main structure of the point cloud [42]. For arbitrary
elements, all field points of each point are searched and the distance from
each point to its neighbor calculated. The mean $\mu$ and the standard
deviation $\sigma$ of the distances of all points from their neighbors are
calculated, and any point whose mean value is outside of the standard interval
will be considered as noise.

$\displaystyle\left[\mu-\alpha_{v}\times\sigma,\
\mu+\alpha_{v}\times\sigma\right]$ (3)

where $\alpha_{v}$<3 determines the width of the standard interval.

The filtered point cloud is shown in Fig. 2 (a).

Figure 2: After color and statistical filter.

####  III-C2 Stitching

According to the hand-eye calibration theory, the orientation of each group of
point clouds (three specific angles) in the coordinate system of the robotic
arm base can be obtained, which provides a good initial orientation for the
following point cloud stitching.

Suppose the number of point cloud points obtained by the camera at specific
angles for $n_{1}$, $n_{2}$ and $n_{3}$ respectively. The coordinates of the k
th point relative to the base coordinate system are ${{}^{R}}P_{A}^{(k)}$,
${{}^{R}}P_{B}^{(k)}$ and ${{}^{R}}P_{C}^{(k)}$ (A, B, C only represent
different positions of the camera) those to the camera coordinate system are
${{}^{CA}}P^{(k)}$, ${{}^{CB}}P^{(k)}$ and ${{}^{CC}}P^{(k)}$ (CA is the
representation of the camera coordinate system at position A, CB, CC are
similar), respectively. Substitute all the values of k to obtain the points
relative to the base coordinate system:

$\left[\begin{array}[]{c}{{}^{R}}P_{A}^{(1)}\\\ {{}^{R}}P_{A}^{(2)}\\\
\vdots\\\ {{}^{R}}P_{A}^{(n_{1})}\end{array}\right]=\ ^{\
R}_{CA}T\left[\begin{array}[]{c}{}^{CA}P^{(1)}\\\ {}^{CA}P^{(2)}\\\ \vdots\\\
{}^{CA}P^{(n_{1})}\end{array}\right]$ (4)

$\left[\begin{array}[]{c}{{}^{R}}P_{B}^{(1)}\\\ {{}^{R}}P_{B}^{(2)}\\\
\vdots\\\ {{}^{R}}P_{B}^{(n_{2})}\end{array}\right]=\ ^{\
R}_{CB}T\left[\begin{array}[]{c}{}^{CB}P^{(1)}\\\ {}^{CB}P^{(2)}\\\ \vdots\\\
{}^{CB}P^{(n_{2})}\end{array}\right]$ (5)

$\left[\begin{array}[]{c}{{}^{R}}P_{C}^{(1)}\\\ {{}^{R}}P_{C}^{(2)}\\\
\vdots\\\ {{}^{R}}P_{C}^{(n_{3})}\end{array}\right]=\ ^{\
R}_{CC}T\left[\begin{array}[]{c}{}^{CC}P^{(1)}\\\ {}^{CC}P^{(2)}\\\ \vdots\\\
{}^{CC}P^{(n_{3})}\end{array}\right]$ (6)

In which, ${}_{CA}^{\ R}T$= ${}_{CB}^{\ R}T$=${}_{CC}^{\
R}T$=${}_{C}^{R}T$=${}_{F}^{R}T^{(i)}$$\widehat{{{}_{C}^{F}}T}$, is the
transformation matrix between the camera coordinate system and the base
coordinate system of the robotic arm deduced in Sec. III-B. Thus, the stitched
point cloud $P_{s}$ can be written as:

$P_{s}=\left[\begin{array}[]{c|c|c}\left[\begin{array}[]{c}{{}^{R}}P_{A}^{(1)}\\\
{{}^{R}}P_{A}^{(2)}\\\ \vdots\\\
{{}^{R}}P_{A}^{(n_{1})}\end{array}\right]^{T}\\!&\left[\begin{array}[]{c}{{}^{R}}P_{B}^{(1)}\\\
{{}^{R}}P_{B}^{(2)}\\\ \vdots\\\
{{}^{R}}P_{B}^{(n_{3})}\end{array}\right]^{T}&\left[\begin{array}[]{c}{{}^{R}}P_{C}^{(1)}\\\
{{}^{R}}P_{C}^{(2)}\\\ \vdots\\\
{{}^{R}}P_{C}^{(n_{3})}\end{array}\right]^{T}\end{array}\right]^{T}$ (7)

The effect of point cloud stitching is shown in Fig. 2 (b).

Figure 3: Fcaf3d network structure.

###  III-D Fcaf3d-lychee for lychee picking point detection

[43] proposed 3D target detection network only using point cloud data as
material, which known as Fcaf3d (Fully Convolutional Anchor-Free 3D Object
Detection). It is a network model that greatly improves detection
generalization specifically designed to achieve fast and accurate detection of
indoor 3D objects. This provides technical support for the identification and
localization of lychee picking points. Fcaf3d network algorithm is an improved
algorithm based on [44] proposed GSDN (Generative sparse detection networks).
Among the improvements, Fcaf3d proposed a fully convolutional anchorless box
and a novel OBB (oriented bounding box) parametrization, which improves<EMAIL_ADDRESS>by 11.4% in comparison with GSDN. Thus, the Fcaf3d network algorithm provides
the best method for the identification and localization of lychees picking
points.

There are four main components of the Fcaf3d network model: Input, Backbone,
Neck, and Head. We modified its structure and combined the SE module with
every Res unit in the backbone to form the SE-Res Block. In this paper, we
will call it Fcaf3d-lychee network model. The structure of the Fcaf3d-lychee
model is shown in Fig. 3.

The sparseness of the point cloud and the fact that the observable surface of
the point cloud does not intersect with the center of the instance, in
addition to the problem of high computational cost of 3d convolution, cause
great difficulties in the task of 3D object detection. Instead of normal dense
3d convolution computation, Fcaf3d and GSDN both adopted the computation
proposed by [45] for sparse convolution of point cloud. The calculation is as
in eq. (8).

Let $\textit{x}_{\textbf{u}}^{in}\in\mathbb{R}^{N^{in}}$ be an
$N^{in}$-dimensional input feature vector in a D-dimensional space at
$\mathbf{u}\in\mathbb{R}^{D}$ (a D-dimensional coordinate), and convolution
kernel weights be
$\textbf{W}\in\mathbb{R}^{{K^{D}}\times{N^{in}}\times{N^{out}}}$. Next, break
down the weights into spatial weights with $K^{D}$ matrices of size
${N^{in}}\times{N^{out}}$ as $\textit{W}_{\textbf{i}}$ for
$|\\{\textbf{i}\\}_{\textbf{i}}|=K^{D}$. Then the generalized sparse
convolution is as follows:

$\textbf{x}^{out}_{\textbf{u}}=\
\sum_{i\in\mathbb{N}^{D}(\textit{u},\mathbb{C}^{in})}W_{\textbf{i}}\textbf{x}_{\textbf{u}+\textbf{i}}^{in}\
\text{for}\ \textbf{u}\in\mathbb{C}^{out}$ (8)

where $\mathbb{N}^{D}$ is a set of offsets that define the shape of a kernel
and $\mathbb{N}^{D}(\textbf{u},\mathbb{C}^{in})=\\{\textbf{i}|\textbf{u+i}\in
C^{in},\textbf{i}\in\mathbb{N}^{D}\\}$ as the set of offsets from the current
center, u, that exist in $\mathbb{C}^{in}$. $\mathbb{C}^{in}$ and
$\mathbb{C}^{out}$ are predefined input and output coordinates of sparse
tensors.

The sparse tensor is different from previous ones in that it is batch-indexed,
i.e., new sparse tensor coordinates are given an index number to distinguish
points from different batches occupying the same coordinates, and it is in COO
format to enable efficient querying of neighbor coordinates. Fcaf3d’s Backbone
adopts a variant form of the same high-dimensional residual structure used in
GSDN, which consists of a series of step-by-step convolutions as well as
residual blocks, each of which reduces the resolution of the target space and
exponentially increases the size of the receptive field. In the Neck, a
transposed convolution is used to improve the features of the large low-
resolution receptive field generated by the previous residual block, and jump-
joining with sparse pruning is utilized to control the sparse tensor of the
cubic growth. The Head, on the other hand, consists of parallel sparse 3d
convolutional layers that output the classification probability, bounding box
regression parameters, and bounding box centrality of the layer, respectively.
On this basis, in order to improve the accuracy and generalization, Danila
Rukhovich et al. proposed the anchor-box-free proposal and the fuzzified
bounding box parameters using the idea of Möbius rings, which greatly reduces
the number of hyper-parameters during training, and the generalization and
accuracy are also dramatically improved.

However, the result of Fcaf3d is undesirable for the identification and
localization of lychee picking points with small targets. Therefore, this
paper proposes the Fcaf3d-lychee model for 3d point cloud detection of picking
points of lychee. For small targets, especially small targets in the point
cloud, the information in the point cloud is inhomogeneous, and the features
of the small targets may be weak with respect to the features in the whole
point cloud, which is very important for the detection of small targets.
Additionally, in the residual structure, the number of input channels per
residual block is too large, while a small number of channels is better for
retaining some special information about small targets in deep features. In
conclusion, the following improvements are made in this paper for the
identification and localization of lychee picking points: (1) Se-block:
Combined the se-attention mechanism module with each res unit in the backbone,
resulting in the SE-Residual module; (2) Hyperparameters: For the four
residual blocks of the resnet, the original input channels are (64, 128, 256,
512), and we get the optimal number of channels for the experiments to be (64,
128, 256, 256), which not only improves the accuracy of the model but also can
reduce the occupation of the video memory greatly.

Figure 4: Structure of SE module.

###  III-E Squeeze-and-Excitation (SE) module

Soft attention mechanism generally includes channel attention mechanism and
spatial attention mechanism. The spatial attention module focuses on the
spatial location information of the data, while the channel attention module
pays more attention to the channel information of the features. SENet belongs
to channel attention, which weights the channels of different features to
establish the relationship between the channels of the features, and the
introduction of SENet will not change the original spatial dimensions. This is
a good choice for sparse and non-uniform data such as point cloud, it just
automatically obtains the weights of the feature channels and then uses the
obtained weights to filter the key features.

The SE module consists of two stages: Squeeze and Excitation. In the Squeeze
stage, the original feature map is dimensionally reduced through global
average pooling, resulting in a single value for each channel. This value
represents the global importance of each channel. In the Excitation stage, a
fully connected layer is used with an input dimension equal to the number of
channels and an output dimension significantly smaller than the number of
channels. The output is then processed through the ReLU activation function to
learn the weights for each channel. These weights are employed to reweight the
features of each channel [46]. The Squeeze-and-Excitation attention mechanism
was expressed as:

$\begin{split}\textbf{$M_{c}$}(F)&=F\cdot\sigma(L_{2}(ReLu(L_{1}(AvgPool(F)))))\\\
&=F\cdot\sigma(L_{2}(ReLu(L_{1}(F_{avg}))))\end{split}$ (9)

where $M_{c}$ is the attention map, the input feature map is F, $\sigma$
represents the sigmoid functions, ReLu is the Non-linear, $L_{1}$ and $L_{2}$
are the full connection operation. $F_{avg}$ is the average pooling operations
on the feature map. The SE attention module is shown in Fig. 4.

###  III-F Lychee Picking Point Dataset

####  III-F1 Data collection

This study was conducted from June 23 to 24, 2023, at lychee Expo Park, East
Huanshi Road, Jiangpu Street, Conghua District, Guangzhou City, Guangdong
Province, China (Longitude $23^{\circ}59^{\prime}$E, Latitude
$113^{\circ}62^{\prime}$N). A total of 850 groups (three point cloud per
group, in preparation for subsequent stitching) of point cloud data were
collected from the dwarf lychee species ”Huai Zhi” using the depth camera
Azure Kinetic DK. These lychee point cloud data will be used to build a lychee
identification and detection dataset. Some invalid point cloud data with
severe occlusion overlap and extremely dark point cloud data were discarded,
leaving 800 point cloud data obtained from different shooting distances as
well as different lighting conditions.

####  III-F2 Voxel Downsampling

As shown in Fig. 5, in this paper, a series of data enhancement operations are
carried out on the original dataset, and data preprocessing is a very
important process, especially for sparse data such as point cloud, the effect
of preprocessing has a great impact on the subsequent production of datasets
and model training, and the preprocessing of the dataset prior to the training
of the model will be described in this subsection. Preprocessing includes but
is not limited to, overall rotation, translation, scaling, and flipping of the
data, but the most important pre-processing steps are multi-view point cloud
alignment and voxel downsampling. These preprocessing operations are very
important in the training of the model, especially for the special data of the
point cloud, these operations improve the robustness of the data, greatly
improve the training accuracy of the model, minimize the overfitting
phenomenon of the model, and improve the generalization ability of the model
training.

Figure 5: Point cloud augmentation.

Although the above pre-processing does not reduce the training difficulty of
the model, the voxel representation is used to further compress the point
cloud data without destroying the original information. In this study, point
cloud downsampling is used in open3d point cloud filtering. The algorithm is
to create a uniformly downsampled point cloud from the input point cloud using
a regular voxel grid in a two-part operation, where first the input point
cloud is voxelised, followed by the centre of mass of the point cloud within
the voxel de-voxelisation of all non-empty voxels as the location of the point
of the final modified voxel. In this way, we further reduce the number of
point clouds, greatly reducing the performance requirements, and also line to
improve the training speed and accuracy, but also improve the final
generalisation of the model. For the choice of voxel size, this study set the
voxel size to 0.01m after many trials and experiments, from the initial two
hundred thousand or so points, reduced to less than one hundred thousand
points, and more importantly, this is the training of the model has a very
obvious increase in the accuracy of training, but also greatly improved the
generalisation of the model. Specific training results are given in the
experimental part of section IV.

##  IV Experiment and Discussion

###  IV-A Dataset and Experimental Setup

In this paper, the specific configuration of the experimental setup is
organized in Table I. The Fcaf3d-litchi model was trained on the mmdetection3d
framework. The specific steps of training are as follows: as in Table II, the
dataset is divided into a training set and a validation set according to the
ratio of 8:2, thus the model is trained on 800 point cloud data divided into
small batches, and the network model adopts the Adamw optimizer, with the
initial learning rate set to 0.001 and the weight decay set to 0.0001. The
same training parameters and platform were also used to train Votenet, Tr3d,
and the initial Fcaf3d models for comparison.

TABLE I: Model training environment settings.

Parameter | Configuration
---|---
Operating system | Ubuntu 20.04
Deep learning framework | 1.8.0+cu113
Programming language | Python 3.8
CUDA | 11.3
CPU | Intel(R) Core(TM) i5-9400F<EMAIL_ADDRESS>GHz
GPU | NVIDIA GeForce RTX 2070 super

TABLE II: Lychee picking points dataset split used for training. | Point cloud | Picking point ground truth
---|---|---
Train | 640 | 1964
Validation | 160 | 516

###  IV-B Metrics

For the evaluation of model training, this study will use, but not be limited
to, metrics such as Precision, Recall, and $F_{1}$ score to assess the
performance of the trained model on the test dataset. The use of the $F_{1}$
score is motivated by its combination of Recall (R) and Precision (P), making
it a comprehensive metric for performance evaluation. Therefore, many previous
studies on fruit detection have widely adopted the $F_{1}$ score as one of the
key evaluation metrics. Accuracy is an important reference parameter to
determine whether a model is performing well or not. TurePositive(TP) is the
number of lychee picking points identified by the model, FalsePositive(FP) is
the number of lychee picking points over-identified by the model,
FalseNegative(FN) is the number of lychee picking points under-identified by
the model.The expressions for Precision, Recall, $F_{1}$ score and Accuracy
are as follows:

$\displaystyle Precision=\frac{TP}{TP+FP}$ (10) $\displaystyle
Recall=\frac{TP}{TP+FN}$ (11) $\displaystyle F_{1}=\frac{2\times
Precision\times Recall}{Precision+Recall}$ (12) $\displaystyle
Accuracy=\frac{TP+TN}{TP+FN+FP+TN}$ (13)

###  IV-C Comparsion of model training results

To verify the recognition and detection performance of the improved
Fcaf3d-lychee model, this study compares it with three classic 3D detection
neural networks: Fcaf3d, Votenet, and Tr3d. Additionally, Fcaf3d-lychee with
different voxel sizes for preprocessing is also compared. Table III presents
the results of these models and the improved Fcaf3d-lychee with varying voxel
sizes for preprocessing.

Figure 6: $F_{1}$ score achieved by four models Fcaf3d, Votenet, Tr3d, and Fcaf3d-lychee model after 300 epochs training. TABLE III: Test results of the 4 models. Model name | Precision(%) | Recall(%) |  $F_{1}$ score(%)
---|---|---|---
Fcaf3d | 60.03 | 66.67 | 63.18
Votenet | 53.31 | 61.11 | 56.94
Tr3d | 63.12 | 74.44 | 68.31
Fcaf3d-lychee | 83.39 | 94.44 | 88.57

As shown in Table III, using the same dataset and voxel preprocessing sizes,
the results of the proposed Fcaf3d-lychee model are noticeable compared to
other detection models. In comparison to Fcaf3d, the Fcaf3d-lychee model
exhibits a 23.36% improvement in Precision, a 27.77% increase in Recall, and a
25.36% rise in $F_{1}$ score. When compared to Votenet, Fcaf3d-lychee shows a
30.08% improvement in Precision, a 33.33% increase in Recall, and a 31.63%
rise in $F_{1}$ score. Compared to Tr3d, the Fcaf3d-lychee model demonstrates
a 20.27% increase in Precision, a 20.00% rise in Recall, and a 20.26%
improvement in $F_{1}$ score. The $F_{1}$ score curves for each model are
depicted in Fig. 6, showing significant improvement in the $F_{1}$ score curve
of the proposed model compared to other models.

###  IV-D Ablation Study

To further validate the improvements proposed in this paper for the
Fcaf3d-lychee model, this section conducts ablation experiments to examine the
performance variations and provide a comprehensive understanding of the
contributions and impact of each improvement on the model’s performance. We
experiment with different voxel sizes, varying input channel configurations in
residual blocks, and comparing models with and without SE attention modules.
The results of the ablation experiments are summarized in Table IV.

####  IV-D1 Voxel Size Preprocessing

The experiments involve using different voxel sizes for preprocessing. The
results indicate that, for the dataset used in this paper, the optimal voxel
size for preprocessing is a voxel cubic grid of 0.01m. Larger voxel sizes lead
to a significant loss of information for small objects in the point cloud,
resulting in a rapid decline in detection performance. On the other hand,
excessively small preprocessing voxel sizes make the overall structure overly
complex, hindering effective detection.

####  IV-D2 SE attention

With optimized voxel preprocessing, experiments are conducted with and without
attention mechanisms. The results show a significant improvement in detection
performance with the introduction of the SE attention mechanism. The channel
attention mechanism with a compression ratio of 16 enhances the detection
performance of the proposed Fcaf3d-lychee model.

####  IV-D3 Res in-channel

After incorporating the SE attention mechanism, modifications are made to the
input channels in residual blocks. The results from the summarized table
demonstrate that adjusting the channel numbers to (64, 128, 256, 256)
significantly improves the model’s performance. The model achieves its best
performance with an accuracy of 77.84% and a recall rate of 94.44%. Optimizing
the channel numbers effectively enhances the model’s performance while
maintaining stability and robustness.

TABLE IV: Ablation Study on Voxel Size, SE-Attention, and Residual Block Input Channels in Fcaf3d-lychee. The better options are marked in bold.Values represent the percentage (%) of Precision, Recall, and $F_{1}$ Score. Ablating parameter | Voxel size(m) | Precision | Recall |  $F_{1}$ score
---|---|---|---|---
Voxel size(m) | NO | 60.03 | 66.67 | 63.18
0.01 | 70.65 | 88.89 | 79.16
SE attention | NO | 70.65 | 88.89 | 79.16
YES | 75.52 | 94.44 | 83.92
Res_inchannel | (64, 128, 256, 512) | 75.52 | 94.44 | 83.92
(64, 128, 256, 256) | 83.39 | 94.44 | 88.57

###  IV-E Visual Sensing in Orchards

Figure 7: Lychee picking robot equipped with Fcaf3d-lychee model in (a)
Guangzhou lychee Expo Park, (b) The success rate of picking point localization
is evaluated by performing cutting action by the end-effector of the picking
robot.

In this section, we show the results of 3D inspection in a lychee garden. This
experiment was conducted on July 2, 2023, from 1-3 pm. The orchard environment
road surface is not very rugged. The data was collected in the afternoon under
strong sunlight. In our previous work, [47] proposed a new semantic mapping
and navigation framework for autonomous navigation of orchard robots. The
picking robot moved on the aisle, and the Fcaf3d-lychee model acquired the
fruit of the target lychee in the real-time point cloud stream as the robot
performed further precise positioning of the picking point location (As shown
in Fig. 7). When the position of the target lychee is obtained, the robotic
arm will be a depth camera suitable depth for three directions of point cloud
acquisition and point cloud stitching in Sec. 3.3, and then the detection
results will be obtained in the Fcaf3d-lychee model.

TABLE V: Performance with different occlusion situations by Fcaf3d-lychee. Occlusion situation of picking points | TP | FP | Accuracy
---|---|---|---
Unobstructed | 232 | 17 | 0.932
Slight obstructed | 103 | 22 | 0.824
Severely obstructed | 75 | 23 | 0.765

To determine the detection effect of the model in real environments, as well
as the accuracy and robustness of the detection for different difficulties,
this experiment was conducted in the case of no occlusion, average occlusion,
and severe occlusion, respectively. For example, Fig. 8(a) shows the case of
no occlusion, Fig. 8(b) shows the case of general occlusion, and Fig. 8(c) and
(d) shows the case of this gear in the eye. The detection performance of the
three cases is shown in Table V, it can be seen that the accuracy of the
picking point in the case of no occlusion is the highest, up to 0.932, and the
highest accuracy of the position, in Fig. 8(c) and (d) this more serious case
of occlusion, the picking point of the gesture will have a greater error, and
there will be interference to make the picking point of the false detection of
the situation occurs, but the position of the point of the picking point is
still very accurate, with an accuracy rate of 0.765. The reason is that under
the occluded eye, the point cloud splicing is not very high quality for the
picking point, which leads to the point cloud distortion at the picking point
or the perception information is less, which makes it difficult for the model
to mistakenly think that it is a picking point, resulting in a low accuracy
rate.

Figure 8: Detection results in the lychee orchard by using the Fcaf3d-lychee.
(a) The occlusion is not serious or there is no occlusion, (b) The situation
of general occlusion, (c) and (d) are the detection results under severe
occlusion.

###  IV-F Discussion

In this paper, a point cloud detection network based on deep learning is
proposed to detect the pose of the target point in color point cloud data
obtained by the depth camera sensor. The point cloud data of the same target
at different three angles (frontal, two horizontal surfaces at 45-degree
angles to the frontal) are stitched using the stitching algorithm, and a large
number of points (200k-300k) after stitched are pre-processed by voxel down-
sampling for efficient target detection of targets with different occlusion
situations. The network detects the pose of the target point on the pre-
processed point cloud data and guides the picking action of the robot arm.

Figure 9: 3D localization errors of the lychee picking point. Each data point
represented the error of that localization in three directions. The larger
cube volume represented a larger localization error.

In the experiments, we analysed the performance impact of different
downsampling sizes, attention mechanisms and the number of input channels of
the residual block on the performance of the proposed Fcaf3d-lychee model.
Experimental results show that the downsampling voxel size has a large impact
on network performance. A voxel size that is too small (less than 0.01m) has
no effect on improving network performance and is not beneficial for model
resource usage. If the voxel size is too large (greater than 0.01m), the
network performance may even be severely degraded due to unclear features. The
experiment also verified the effect of the SE attention mechanism and the
number of input channels of the residual block on the performance of the
proposed model. Overall, with an appropriate choice of sampling voxel size,
multi-viewpoint cloud stitching, attention mechanism fusion and number of
input channels, the proposed Fcaf3d-lychee model achieves a $F_{1}$ score of
88.57% on the test dataset. The target detection accuracy under different
occlusion conditions in the real orchard environment reached 0.932, 0.824 and
0.765, respectively. To more simply express the 3D localisation error of
Fcaf3d-lychee on the lychee picking point. We performed 100 localisation
accuracy experiments in different scenarios. Then, the 3D error graph was
plotted (Fig. 9), which shows that the localisation errors in three directions
are within the range of ±1.5 cm, reflecting the high localisation accuracy of
the fcaf3d-lychee proposed in this paper.

##  V Conclusion

This paper presents a novel method based on an enhanced 3D point cloud object
detection network to accurately identify lychee-picking points in unstructured
natural environments. The aim is to determine optimal picking points to guide
picking robots to harvest lychee at suitable angles. Based on the Fcaf3d
model, the research focuses on training lychee-picking point detection models
to achieve high accuracy, robustness and generalisation. The presented new
detection network Fcaf3d-lychee, which incorporates attention and dynamic
sparse pruning, shows excellent performance in various scenarios with an
average accuracy (AP) of 83.39% and a $F_{1}$ value of 88.57. Compared to
other 3D detection models such as Tr3d and Votenet, the proposed network
achieves significantly higher accuracy. Unlike traditional 2D detection
networks such as YOLO, the network proposed in this paper directly obtains the
3D coordinates of lychee picking points, achieving effects that were
previously unattainable. Experimental results in real orchards show that the
Fcaf3d-lychee model achieves robust detection performance under different
shading conditions, with accuracies of 0.932, 0.824, and 0.765 for no shading,
general shading, and heavy shading, respectively. Furthermore, the
localisation accuracy of the Fcaf3d-lychee remains within ±1.5 cm in all
directions, which is within the acceptable range of picking error. The
proposed Fcaf3d-lychee target detection network overcomes the limitations of
previous methods, lays the foundation for the accurate calculation of picking
points, and provides essential technical support for the practical application
of lychee picking robots.

## References

  * [1] H. Zhou, X. Wang, W. Au, H. Kang, and C. Chen, “Intelligent robots for fruit harvesting: recent developments and future challenges,” _PRECISION AGRICULTURE_ , vol. 23, no. 5, pp. 1856–1907, OCT 2022.
  * [2] Z. Zhong, J. Xiong, Z. Zheng, B. Liu, S. Liao, Z. Huo, and Z. Yang, “A method for litchi picking points calculation in natural environment based on main fruit bearing branch detection,” _Computers and Electronics in Agriculture_ , vol. 189, p. 106398, 2021.
  * [3] X. Wang, H. Kang, H. Zhou, W. Au, and C. Chen, “Geometry-aware fruit grasping estimation for robotic harvesting in apple orchards,” _Computers and Electronics in Agriculture_ , vol. 193, p. 106716, 2022.
  * [4] C. Bac, J. Hemming, and E. van Henten, “Stem localization of sweet-pepper plants using the support wire as a visual cue,” _Computers and Electronics in Agriculture_ , vol. 105, pp. 111–120, 2014.
  * [5] G. Coll-Ribes, I. J. Torres-Rodríguez, A. Grau, E. Guerra, and A. Sanfeliu, “Accurate detection and depth estimation of table grapes and peduncles for robot harvesting, combining monocular depth estimation and cnn methods,” _Computers and Electronics in Agriculture_ , vol. 215, p. 108362, 2023.
  * [6] J. Xiong, R. Lin, Z. Liu, Z. He, L. Tang, Z. Yang, and X. Zou, “The recognition of litchi clusters and the calculation of picking point in a nocturnal natural environment,” _BIOSYSTEMS ENGINEERING_ , vol. 166, pp. 44–57, FEB 2018.
  * [7] L. Luo, Y. Tang, Q. Lu, X. Chen, P. Zhang, and X. Zou, “A vision methodology for harvesting robot to detect cutting points on peduncles of double overlapping grape clusters in a vineyard,” _COMPUTERS IN INDUSTRY_ , vol. 99, pp. 130–139, AUG 2018.
  * [8] X. Liu, D. Zhao, W. Jia, W. Ji, and Y. Sun, “A detection method for apple fruits based on color and shape features,” _IEEE Access_ , vol. 7, pp. 67 923–67 933, 2019.
  * [9] Y. Bai, S. Mao, J. Zhou, and B. Zhang, “Clustered tomato detection and picking point location using machine learning-aided image analysis for automatic robotic harvesting,” _PRECISION AGRICULTURE_ , vol. 24, no. 2, pp. 727–743, APR 2023.
  * [10] L. Huang, S. X. Yang, and D. He, “Abscission point extraction for ripe tomato harvesting robots,” _Intelligent Automation & Soft Computing_, vol. 18, no. 6, pp. 751–763, 2012.
  * [11] G. Lin, Y. Tang, X. Zou, J. Cheng, and J. Xiong, “Fruit detection in natural environment using partial shape matching and probabilistic hough transform,” _Precision Agriculture_ , vol. 21, pp. 160–177, 2020.
  * [12] C. Chen, J. Lu, M. Zhou, J. Yi, M. Liao, and Z. Gao, “A yolov3-based computer vision system for identification of tea buds and the picking point,” _COMPUTERS AND ELECTRONICS IN AGRICULTURE_ , vol. 198, JUL 2022.
  * [13] H. Wang, Y. Lin, X. Xu, Z. Chen, Z. Wu, and Y. Tang, “A study on long-close distance coordination control strategy for litchi picking,” _AGRONOMY-BASEL_ , vol. 12, no. 7, JUL 2022.
  * [14] Z. Xiong, L. Wang, Y. Zhao, and Y. Lan, “Precision detection of dense litchi fruit in uav images based on improved yolov5 model,” _Remote Sensing_ , vol. 15, no. 16, 2023.
  * [15] C. Li, J. Lin, Z. Li, C. Mai, R. Jiang, and J. Li, “An efficient detection method for litchi fruits in a natural environment based on improved yolov7-litchi,” _Computers and Electronics in Agriculture_ , vol. 217, p. 108605, 2024.
  * [16] H. Peng, C. Xue, Y. Shao, K. Chen, H. Liu, J. Xiong, H. Chen, Z. Gao, and Z. Yang, “Litchi detection in the field using an improved yolov3 model,” _International Journal of Agricultural and Biological Engineering_ , vol. 15, no. 2, pp. 211–220, 2022.
  * [17] J. Li, Y. Tang, X. Zou, G. Lin, and H. Wang, “Detection of fruit-bearing branches and localization of litchi clusters for vision-based harvesting robots,” _IEEE Access_ , vol. 8, pp. 117 746–117 758, 2020.
  * [18] H. Kang, X. Wang, and C. Chen, “Accurate fruit localisation using high resolution lidar-camera fusion and instance segmentation,” _Computers and Electronics in Agriculture_ , vol. 203, p. 107450, 2022.
  * [19] X. Du, Z. Meng, Z. Ma, W. Lu, and H. Cheng, “Tomato 3d pose detection algorithm based on keypoint detection and point cloud processing,” _Computers and Electronics in Agriculture_ , vol. 212, p. 108056, 2023.
  * [20] T. Li, Q. Feng, Q. Qiu, F. Xie, and C. Zhao, “Occluded apple fruit detection and localization with a frustum-based point-cloud-processing approach for robotic harvesting,” _REMOTE SENSING_ , vol. 14, no. 3, FEB 2022.
  * [21] L. Zeng, J. Feng, and L. He, “Semantic segmentation of sparse 3d point cloud based on geometrical features for trellis-structured apple orchard,” _Biosystems Engineering_ , vol. 196, pp. 46–55, 2020.
  * [22] T. Yu, C. Hu, Y. Xie, J. Liu, and P. Li, “Mature pomegranate fruit detection and location combining improved f-pointnet with 3d point cloud clustering in orchard,” _COMPUTERS AND ELECTRONICS IN AGRICULTURE_ , vol. 200, SEP 2022.
  * [23] J. Li, M. Wu, and H. Li, “3d reconstruction and volume estimation of jujube using consumer-grade rgb-depth sensor,” _IEEE ACCESS_ , vol. 11, pp. 61 502–61 512, 2023.
  * [24] H. Kang and X. Wang, “Semantic segmentation of fruits on multi-sensor fused data in natural orchards,” _COMPUTERS AND ELECTRONICS IN AGRICULTURE_ , vol. 204, JAN 2023.
  * [25] T. Yoshida, T. Fukao, and T. Hasegawa, “Cutting point detection using a robot with point clouds for tomato harvesting,” _Journal of Robotics and Mechatronics_ , vol. 32, no. 2, pp. 437–444, 2020.
  * [26] G.-Z. Yang, J. Bellingham, P. E. Dupont, P. Fischer, L. Floridi, R. Full, N. Jacobstein, V. Kumar, M. McNutt, and R. Merrifield, “The grand challenges of science robotics,” _Science Robotics_ , vol. 3, 2018.
  * [27] Y. Li, Q. Feng, T. Li, F. Xie, C. Liu, and Z. Xiong, “Advance of target visual information acquisition technology for fresh fruit robotic harvesting: A review,” _AGRONOMY-BASEL_ , vol. 12, no. 6, JUN 2022.
  * [28] Y. Tang, M. Chen, C. Wang, L. Luo, J. Li, G. Lian, and X. Zou, “Recognition and localization methods for vision-based fruit picking robots: A review,” _frontiers in piant science_ , vol. 11, 2020.
  * [29] Y. Zhao, L. Gong, Y. Huang, and C. Liu, “A review of key techniques of vision-based control for harvesting robot,” _Computers and Electronics in Agriculture_ , vol. 127, pp. 311–323, 2016.
  * [30] L. Li, Z. Ouyang, W. Tian, and W. Sun, “Location and model reconstruction algorithm for overlapped and sheltered spherical fruits based on geometry,” _International Journal of Advanced Robotic Systems_ , vol. 19, p. 172988142110567, 2022.
  * [31] Z.-L. He, J.-T. Xiong, R. Lin, X. Zou, L.-Y. Tang, Z.-G. Yang, Z. Liu, and G. Song, “A method of green litchi recognition in natural environment based on improved lda classifier,” _Computers and Electronics in Agriculture_ , vol. 140, pp. 159–167, 2017.
  * [32] M. Afonso, H. Fonteijn, F. S. Fiorentin, D. Lensink, M. Mooij, N. Faber, G. Polder, and R. Wehrens, “Tomato fruit detection and counting in greenhouses using deep learning,” _Frontiers in Plant Science_ , vol. 11, 2020.
  * [33] C. Liang, J. Xiong, Z. Zheng, Z. Zhong, Z. Li, S. Chen, and Z. Yang, “A visual detection method for nighttime litchi fruits and fruiting stems,” _Computers and Electronics in Agriculture_ , vol. 169, p. 105192, 2020.
  * [34] B. Yan, P. Fan, X. Lei, Z. Liu, and F. Yang, “A real-time apple targets detection method for picking robot based on improved yolov5,” _REMOTE SENSING_ , vol. 13, no. 9, MAY 2021.
  * [35] Y. Jin, J. Liu, J. Wang, Z. Xu, and Y. Yuan, “Far-near combined positioning of picking-point based on depth data features for horizontal-trellis cultivated grape,” _Computers and Electronics in Agriculture_ , vol. 194, p. 106791, 2022.
  * [36] C. Zhang, K. Zhang, L. Ge, K. Zou, S. Wang, J. Zhang, and W. Li, “A method for organs classification and fruit counting on pomegranate trees based on multi-features fusion and support vector machine by 3d point cloud,” _Scientia Horticulturae_ , vol. 278, p. 109791, 2021.
  * [37] Z. Zhang, “A flexible new technique for camera calibration,” _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , vol. 22, no. 11, pp. 1330–1334, 2000.
  * [38] M. Chen, Y. Tang, X. Zou, Z. Huang, H. Zhou, and S. Chen, “3d global mapping of large-scale unstructured orchard integrating eye-in-hand stereo vision and slam,” _Computers and Electronics in Agriculture_ , vol. 187, p. 106237, 2021.
  * [39] Y. Tang, L. Li, C. Wang, M. Chen, W. Feng, X. Zou, and K. Huang, “Real-time detection of surface deformation and strain in recycled aggregate concrete-filled steel tubular columns via four-ocular vision,” _Robotics and Computer-Integrated Manufacturing_ , vol. 59, pp. 36–46, 2019.
  * [40] H. Kang, Y. Zang, X. Wang, and Y. Chen, “Uncertainty-driven spiral trajectory for robotic peg-in-hole assembly,” _IEEE Robotics and Automation Letters_ , vol. 7, no. 3, pp. 6661–6668, 2022.
  * [41] M. Chunyan, Z. Lihua, S. Hong, and Y. Wei, “Research on 3d reconstruction of fruit tree and fruit recognition and location method based on rgb-d camera,” _Nongye Jixie Xuebao/Transactions of the Chinese Society of Agricultural Machinery_ , 2015.
  * [42] M. Chen, Y. Tang, X. Zou, K. Huang, L. Li, and Y. He, “High-accuracy multi-camera reconstruction enhanced by adaptive point cloud correction algorithm,” _Optics and Lasers in Engineering_ , vol. 122, pp. 170–183, 2019.
  * [43] D. Rukhovich, A. Vorontsova, and A. Konushin, “Fcaf3d: Fully convolutional anchor-free 3d object detection,” _Computer Vision – ECCV 2022_ , pp. 477–493, 2022.
  * [44] J. Gwak, C. Choy, and S. Savarese, “Generative sparse detection networks for 3d single-shot object detection,” _Computer Vision – ECCV 2020_ , pp. 297–313, 2020.
  * [45] C. Choy, J. Gwak, and S. Savarese, “4d spatio-temporal convnets: Minkowski convolutional neural networks,” _2019 IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR)_ , 2019.
  * [46] J. Hu, L. Shen, and G. Sun, “Squeeze-and-excitation networks,” _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ , June 2018.
  * [47] Y. Pan, K. Hu, H. Cao, H. Kang, and X. Wang, “A novel perception and semantic mapping method for robot autonomy in orchards,” _Computers and Electronics in Agriculture_ , vol. 219, p. 108769, 2024.


#  Gluon and valence quark distributions for the pion and kaon in nuclear
matter

Parada T.P. Hutauruk<EMAIL_ADDRESS>Department of Physics, Pukyong
National University (PKNU), Busan 48513, Korea  Seung-il Nam<EMAIL_ADDRESS>Department of Physics, Pukyong National University (PKNU), Busan 48513, Korea
Center for Extreme Nuclear Matters (CENuM), Korea University, Seoul 02841,
Korea

###### Abstract

In this paper we study the gluon and valence quark distributions in the pion
and kaon in nuclear medium for various nuclear densities as well as in vacuum
within the Nambu–Jona-Lasinio (NJL) model with the help of the proper-time
regularization scheme which simulates a confinement of QCD. The nuclear medium
effect is also determined in the same model for the symmetric nuclear matter.
We then analyze the gluon and valence quark distributions for the bound pion
and kaon in symmetric nuclear matter as well as those in vacuum. We find that
the valence quark and gluon distributions in vacuum have relatively good
agreements with the experimental data, the lattice QCD simulations, and the
JAM Monte-Carlo (MC) global fit QCD analysis. Evolving to the higher
factorization scale $Q=$ 4 GeV, the in-medium gluon and valence-quark
distributions of the pion for various nuclear densities are turned out to be
almost unchanged in comparison to the vacuum cases. On the contrary, for the
kaon, they increase significantly with respect to the densities. Finally, we
find that the vacuum gluon distribution for the kaon is smaller than that for
the pion, which is consistent with other theoretical predictions. This feature
holds for the in-medium gluon distribution in the nuclear density up to the
saturation density.

##  I Introduction

Parton distribution functions (PDFs) play important roles in describing the
non-perturbative aspects of Quantum ChromoDynamics (QCD) for the internal
structure of the hadron bound states [1]. Also, they are very crucial
quantities for computing the cross sections for the high-energy hadron-hadron
and neutrino-hadron interactions. Recently the valence quark and gluon
distributions of the pseudo-scalar (PS) mesons, i.e., kaon and pion become a
more attractive subject in hadronic physics community [2, 3, 4, 5, 6, 7, 8, 9,
10]. This is due to they can provide us a better understanding for the PS
mesons’internal structures and the dynamics of quarks and gluons. In addition
to these theoretical reasons, it is also triggered by the future experiments,
i.e, the Electron-Ion Collider (EIC) [6, 11] and Electron-Ion Collider in
China (EicC) [12] as well as the COMPASS++/AMBER new QCD facility at CERN SPS
[13], which will be expected to provide more precise data with a wide range of
the kinematic coverage for the gluon and quark distributions for the mesons.

It is widely known that pions and kaons are the Nambu-Goldstone bosons (NGBs)
that emerged as a consequence of dynamical chiral symmetry breaking
(D$\chi$SB). A deeper understanding of the internal structure of the mesons
leads us to more profound insights of the D$\chi$SB and  vice versa. Recently,
remarkable progresses have been achieved in studying the gluon distributions
of the PS mesons in vacuum. Several theoretical models [4, 8, 10, 14], the
Jefferson Lab Angular Momentum Collaboration (JAM) phenomenology global fit
QCD analysis [5], and lattice QCD simulations [15, 16] have been applied to
study the valence-quark and gluon distribution functions for the mesons in
vacuum. However, besides these efforts, a theoretical understanding of those
quark and gluon distributions still requires more studies, and the situation
is worsened by the scarce of experimental data that brings us to a difficulty
to resolve the current controversies on the valence-quark distributions at
high-$x$, for example, the power law of the distribution functions (DFs) at
$x\rightarrow 1$ [17, 18] and the gluon distributions at low $x$ [18],
respectively. Attempts on studying the valence quark and gluon distribution
functions in nuclear medium are much more limited.

Recent studies on the in-medium modifications of pion and kaon structures have
been reported in the literature [19, 20] for investigating the pion valence-
quark distribution amplitude (VDA) in medium using the light-front (LF) model,
associated with the quark-meson coupling (QMC) model [19] and the in-medium
valence-quark distribution functions (VDFs) of the pion and kaon in the
Nambu–Jona-Lasinio (NJL) model augmented with the QMC model [20]. Similarly,
but for the nucleon case, a very recent study was made to investigate the
gluon distributions of nucleon in vacuum compared with those in the nuclei
using the NJL model [21]. They found the significant effects of the medium
modifications for the unpolarized and polarized gluon distributions of bound
nucleon in nuclear matter. Also, they reported that, in their work, the gluon
distributions of the nucleon at scale $Q^{2}$ were dynamically generated via
next-to-leading order (NLO) Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP)
QCD evolution [22], meaning that it is no gluon dynamics at initial scale.

Inspired by these recent studies, for the first time, we study, in the present
work, the pion and kaon gluon distributions in nuclear medium in the framework
of the NJL model with the help of the proper-time regularization (PTR) scheme,
simulating a confinement of QCD. Aside from focusing on the pion and kaon
gluon distributions in nuclear medium, we also compute the pion and kaon
valence-quark distributions in nuclear medium as well. The NJL model has been
successfully applied for various physics phenomena of low-energy non-
perturbative QCD, namely, the transverse momentum dependent (TMD) [23], the
fragmentation function (FF) [24], the pion and kaon properties in nuclear
medium [25, 26] and the properties of neutron star [27].

In our approach, the nuclear medium effect is also calculated in the NJL
model, which is the same model as used in vacuum one. Thus, we firstly present
our results for the pion and kaon gluon distributions in vacuum compared with
the lattice QCD simulations [15, 16] and the JAM phenomenology global fit QCD
analysis [5]. Note that the gluon distributions for the pion and kaon are
absent at initial model scale as in Ref. [21], since there is no gluon
dynamics in the NJL model. It is absorbed into the $G_{\pi}$ coupling
constant. Hence, in this work, we purely generate the gluon and sea-quark
distributions from the NLO DGLAP QCD evolution through the parton splitting
functions. Next, we compute the gluon and valence quark distributions for the
mesons in nuclear medium for various baryon densities to observe how the in-
medium modifications change the distributions. This study is very helpful to
shed light on the quark-gluon dynamics in the bound pion and kaon in symmetric
nuclear matter (SNM).

This paper is organized as follows. In Sec. II we briefly introduce the
effective Lagrangian for the SU(3) flavor NJL model and the quark properties,
namely, the constituent quark mass and meson-quark coupling constant that
required in the calculation of the vacuum parton distribution functions. We
then present an expression for the twist-2 valence quark distributions of the
pion and kaon in vacuum. Moreover, we describe the SNM of the NJL model (SNM-
NJL) that used for computing the valence quark distributions in SNM. Finally,
we present our formula for the in-medium distributions in the NJL model. In
Sec. III our numerical results are presented and their implications are
discussed. Sec. IV is devoted for a summary.

##  II Nambu–Jona-Lasinio Model

In this section we briefly present the NJL effective Lagrangian and the
properties of the pion and kaon in the model. It maintains the important
features of the non-perturbative QCD, i.e. the spontaneously chiral symmetry
breaking (S$\chi$SB) for instance. The SU(3) flavor NJL Lagrangian is given by

$\displaystyle\mathscr{L}_{\textrm{NJL}}$ $\displaystyle=$
$\displaystyle\bar{q}\left(i\partial\\!\\!\\!/-\hat{m}_{q}\right)q+G_{\pi}[(\bar{q}\mathbf{\lambda}_{a}q)^{2}-(\bar{q}\mathbf{\lambda}_{a}\gamma_{5}q)^{2}]$
(1) $\displaystyle-$ $\displaystyle
G_{\rho}[(\bar{q}\mathbf{\lambda}_{a}\gamma^{\mu})^{2}+(\bar{q}\mathbf{\lambda}_{a}\gamma^{\mu}\gamma_{5}q)^{2}].$

With the quark fields $q$ are defined by $q=(u,d,s)^{T}$,
$\hat{m}_{q}=\textrm{diag}(m_{u},m_{d},m_{s})$ represents the current-quark
mass matrix, and $\mathbf{\lambda}_{a}$ are the Gell-Mann matrices in flavor
space with $\lambda_{0}\equiv\sqrt{\frac{2}{3}}\mathds{1}$. The $G_{\pi}$ and
$G_{\rho}$ are the coupling constants of the four-fermion dimensional with
units of GeV-2. Thus, the standard solution to NJL gap equation is given by

$\displaystyle S_{q}^{-1}(p)$ $\displaystyle=$ $\displaystyle
p\\!\\!\\!/-M_{q}+i\epsilon,$ (2)

where the subscript $q=(u,d,s)$ denotes the quark flavor and the dynamical
quark mass $M_{q}$ in the PTR scheme is given by

$\displaystyle M_{q}$ $\displaystyle=$ $\displaystyle
m_{q}+\frac{3G_{\pi}M_{q}}{\pi^{2}}\int_{\tau_{\textrm{UV}}^{2}}^{\tau_{\textrm{IR}}^{2}}\frac{d\tau}{\tau^{2}}\exp\left(-\tau
M_{q}^{2}\right),$ (3)

where $\tau_{\textrm{IR}}^{2}=1/\Lambda_{\textrm{IR}}^{2}$ and
$\tau_{\textrm{UV}}^{2}=1/\Lambda_{\textrm{UV}}^{2}$ stand for respectively
the infrared (IR) and ultraviolet (UV) integration limits with the value of
$\Lambda_{\textrm{IR}}=$ 0.240 GeV, which determined based on the limit of
$\Lambda_{\textrm{QCD}}$ ($\simeq$ 0.2-0.3 GeV), and $\Lambda_{\textrm{UV}}$
are the infrared and ultraviolet cutoffs, respectively.

In the NJL model, pions and kaons, as bound state of the dressed quark-
antiquark, can be obtained by solving the Bethe-Salpeter equations (BSEs). The
BSE solutions are given by the interaction channel of the two-body amplitude.
For those PS mesons, it has the form

$\displaystyle t_{\pi,K}$ $\displaystyle=$
$\displaystyle\frac{-2iG_{\pi}}{1+2G_{\pi}\Pi_{\pi,K}(p^{2})},$ (4)

where the polarization insertions for the PS mesons are respectively given by

$\displaystyle\Pi_{\pi}(p^{2})$ $\displaystyle=$ $\displaystyle
6i\int\frac{d^{4}k}{(2\pi)^{4}}[\gamma_{5}S_{l}(k)\gamma_{5}S_{l}(k+p)],$ (5)
$\displaystyle\Pi_{K}(p^{2})$ $\displaystyle=$ $\displaystyle
6i\int\frac{d^{4}k}{(2\pi)^{4}}[\gamma_{5}S_{l}(k)\gamma_{5}S_{s}(k+p)].$ (6)

Here, the subscripts of $l=(u,d$) and $s$ are respectively referring to the
light and strange quarks. From the pole of the amplitude $t_{\pi,K}$ in Eq.
(4) we determine the PS-meson masses by solving the pole equations
$1+2G_{\pi}\Pi_{\pi}(p^{2}=m_{\pi}^{2})=$ 0 for the pion and
$1+2G_{\pi}\Pi_{K}(p^{2}=m_{K}^{2})=$ 0 for the kaon. By solving these
equations analytically, the expressions for the PS-meson masses in the PTR
scheme are straightforwardly obtained by

$\displaystyle m_{\pi}^{2}$ $\displaystyle=$
$\displaystyle\frac{m_{l}}{M_{l}}\frac{2}{G_{\pi}\mathscr{I}_{ll}(m_{\pi}^{2})},$
(7) $\displaystyle m_{K}^{2}$ $\displaystyle=$
$\displaystyle\left(\frac{m_{s}}{M_{s}}+\frac{m_{l}}{M_{l}}\frac{1}{G_{\pi}\mathscr{I}_{ls}(m_{K}^{2})}+(M_{s}-M_{l})^{2}\right),$
(8)

where the quantity of $\mathscr{I}_{ls}(p^{2})$ in Eq. (8) is defined for the
kaon by

$\displaystyle\mathscr{I}_{ls}(p^{2})$ $\displaystyle=$
$\displaystyle\frac{3}{\pi^{2}}\int_{0}^{1}dz\int_{\tau_{\textrm{UV}^{2}}}^{\tau_{\textrm{IR}^{2}}}\frac{d\tau}{\tau}$
(9) $\displaystyle\times$
$\displaystyle\exp\left(-\tau(z(z-1)p^{2}+zM_{s}^{2}+(1-z)M_{l}^{2})\right).$

For the pion case, one of the quark flavors is changed as $s\to l$, giving
$\mathscr{I}_{ll}(p^{2})$.

Thus, straightforwardly, the meson-quark coupling constants can be determined
through the residue at pole in the quark-antiquark amplitude $t_{\pi,K}$. In
other word, it can be simply obtained by calculating the first derivative of
the polarization insertion in Eqs. (5), (6) for the corresponding PS meson
with respect to the $p^{2}$. It has a form

$\displaystyle g_{mq\bar{q}}^{-2}$ $\displaystyle=$
$\displaystyle-\frac{\partial\Pi_{m}(p^{2})}{\partial
p^{2}}\Big{|}_{p^{2}=m_{m}^{2}},$ (10)

with the subscript $m=(\pi,K)$. The meson-quark coupling constant of the
$g_{mq\bar{q}}$ has a relation with the wave function renormalization constant
that gives $Z_{m}=g_{mq\bar{q}}^{-2}$.

###  II.1 Vacuum pion and kaon parton distributions

In this section we present the formulas for the valence quark distributions of
the PS mesons in vacuum. The twist-2 quark distributions are simply defined by

$\displaystyle q_{m}(x)$ $\displaystyle=$ $\displaystyle\frac{p^{+}}{2\pi}\int
d\xi^{-}\exp\left(ixp^{+}\xi^{-}\right)\langle
m|\bar{q}(0)\gamma^{+}q(\xi^{-})|m\rangle_{c},$

where $x=\frac{k^{+}}{p^{+}}$ is the Bjorken scaling variable or the
longitudinal momentum of the parton in the PS mesons with $k^{+}$ is the plus-
component of the struck momentum of quark and $p^{+}$ is the plus-component of
the PS meson momentum, $\xi$ is the skewness variable, and the subscript $c$
denotes the connected matrix element.

Figure 1:  (Color online) The relevant diagrams for the PS-mesons valence
quark distributions. The red crossed is an operator insertion
$\gamma^{+}\delta\left(p^{+}x-k^{+}\right)\hat{P}_{q}$, where $\hat{P}_{q}$ is
the projection operator for quarks of flavor $q$.

Following our previous work in Ref. [28], we then calculate the valence quark
distributions based on two Feynman diagrams in Fig. 1. The expressions for the
operator insertion of the light and strange quarks are given as follows

$\displaystyle\gamma^{+}\delta\left(k^{+}-xp^{+}\right)\hat{P}_{l}$
$\displaystyle=$
$\displaystyle\gamma^{+}\delta\left(k^{+}-xp^{+}\right)\frac{1}{2}\left(\frac{2}{3}\mathds{1}\pm\lambda_{3}+\frac{1}{\sqrt{3}}\lambda_{8}\right),$
(12) $\displaystyle\gamma^{+}\delta\left(k^{+}-xp^{+}\right)\hat{P}_{s}$
$\displaystyle=$
$\displaystyle\gamma^{+}\delta\left(k^{+}-xp^{+}\right)\left(\frac{1}{3}\mathds{1}-\frac{1}{\sqrt{3}}\lambda_{8}\right).$
(13)

Considering the relation $\bar{q}(x)=-q(-x)$, the valence-quark distributions
can be defined by

$\displaystyle q_{\pi,K}(x)$ $\displaystyle=$ $\displaystyle
ig_{mq\bar{q}}^{2}\int\frac{d^{4}k}{(2\pi)^{4}}\delta\left(k^{+}-xp^{+}\right)$
$\displaystyle\times$
$\displaystyle\mathrm{Tr}_{c,f,\gamma}[\gamma_{5}\lambda_{a}^{\dagger}S_{l}(k)\gamma^{+}\hat{P}_{u/d}S_{l}(k)\gamma_{5}\lambda_{a}S_{l}(k-p)],$
$\displaystyle\bar{q}_{\pi,K}(x)$ $\displaystyle=$ $\displaystyle-
ig_{mq\bar{q}}^{2}\int\frac{d^{4}k}{(2\pi)^{4}}\delta\left(k^{+}+xp^{+}\right)$
$\displaystyle\times$
$\displaystyle\mathrm{Tr}_{c,f,\gamma}[\gamma_{5}\lambda_{a}S_{l}(k)\gamma^{+}\hat{P}_{\bar{d}/\bar{s}}S_{l}(k)\gamma_{5}\lambda_{a}^{\dagger}S_{s}(k+p)],$

where the meson-quark coupling constant of the $g_{mq\bar{q}}$ that defined in
Eq. (10). The trace runs over the color, flavor, and Lorentz indices. Similar
with Ref. [28], these valence quark distributions are evaluated using the
moment, which is defined by

$\displaystyle\mathscr{A}_{n}$ $\displaystyle=$
$\displaystyle\int_{0}^{1}dxx^{(n-1)}q_{\pi,K}(x),$ (15)

with $n$ is an integer number. We then apply the Ward-Takahashi-like identity
$S(k)\gamma^{+}S(k)=-\partial S(k)/\partial k_{+}$ and perform the Feynman
parameterization.

For general interest, here, we present the explicit expression for the valence
quark distributions of the kaon in the PTR scheme:

$\displaystyle q_{K}(x)$ $\displaystyle=$
$\displaystyle\frac{3g_{Kq\bar{q}^{2}}}{4\pi^{2}}\int d\tau$
$\displaystyle\times$
$\displaystyle\exp\left(-\tau(x(x-1)m_{K}^{2}+xM_{s}^{2}+(1-x)M_{l}^{2})\right)$
$\displaystyle\times$
$\displaystyle\Big{[}\frac{1}{\tau}+x(1-x)\left(m_{K}^{2}-(M_{l}-M_{s})^{2}\right)\Big{]},$

and for the valence antiquark distribution of the kaon is given by

$\displaystyle\bar{q}_{K}(x)$ $\displaystyle=$
$\displaystyle\frac{3g_{Kq\bar{q}^{2}}}{4\pi^{2}}\int d\tau$ (17)
$\displaystyle\times$
$\displaystyle\exp\left(-\tau(x(x-1)m_{K}^{2}+xM_{l}^{2}+(1-x)M_{s}^{2})\right)$
$\displaystyle\times$
$\displaystyle\Big{[}\frac{1}{\tau}+x(1-x)\left(m_{K}^{2}-(M_{l}-M_{s})^{2}\right)\Big{]}.$

Thus, that of the pion is straightforwardly obtained by replacing
$M_{s}\rightarrow M_{l}$ and $g_{Kq\bar{q}}^{2}\rightarrow g_{\pi
q\bar{q}}^{2}$ in Eqs. (II.1), (17), resulting the
$u_{\pi}(x)=\bar{d}_{\pi}(x)$.

The valence quark distributions must satisfy the baryon number and momentum
sum rules. The expression for the baryon number and momentum sum rules for the
PS mesons are respectively given by

$\displaystyle\int_{0}^{1}dx[u_{K}(x)-\bar{u}_{K}(x)]=\int_{0}^{1}dx[\bar{s}_{K}-s_{K}(x)]$
$\displaystyle=$ $\displaystyle 1,$
$\displaystyle\int_{0}^{1}dx\,x[u_{K}(x)+\bar{u}_{K}(x)+s_{K}(x)+\bar{s}_{K}(x)]$
$\displaystyle=$ $\displaystyle 1.$ (18)

Here we again emphasize that, in our NJL model, the PS-meson sea-quark and
gluon distributions are zero at initial scale $Q_{0}^{2}$, since in the NJL
model the gluons are absorbed into the $G_{\pi}$ coupling constant. Hence the
NJL model model has no dynamical gluons as in Ref. [21]. So, they are purely
generated using the NLO DGLAP QCD evolution in present work.

In the DGLAP QCD evolution, the quark distributions for the PS meson are
classified into three types; the gluon distribution, the quark distribution,
and the sea-quark distribution. The quark distribution, which is well known as
the valence quark distribution, is the non-singlet (NS) quark distribution,
which is simply defined by

$\displaystyle q_{\mathrm{NS}}(x)$ $\displaystyle=$ $\displaystyle
q(x)-\bar{q}(x),$ (19)

where $q(x)$ and $\bar{q}(x)$ are respectively the quark and antiquark
distributions. The evolution of the non-singlet quark distribution in the
DGLAP QCD evolution are defined by

$\displaystyle\frac{\partial q_{\textrm{NS}}(x,Q^{2})}{\partial\ln(Q^{2})}$
$\displaystyle=$
$\displaystyle\mathscr{P}_{qq}\left(x,\alpha_{\textrm{s}}(Q^{2})\right)\otimes
q_{\textrm{NS}}(x,Q^{2}),$ (20)

where $\mathscr{P}_{qq}$ denotes the splitting of the q-q (quark-quark)
function. The physics interpretation of $\mathscr{P}_{qq}$ is the probability
for a quark of type $q$ with momentum fraction $z$ to emit the quark and
becomes a quark of new type $q$ with momentum fraction $x$. A convolution
product between the splitting function and the non-singlet quark distribution
has a form

$\displaystyle\mathscr{P}_{qq}\otimes q_{\mathrm{NS}}$ $\displaystyle=$
$\displaystyle\int_{x}^{1}\frac{dz}{x}\mathscr{P}\left(\frac{x}{z}\right)q_{\mathrm{NS}}(z,Q^{2}).$
(21)

The other type of quark distribution is called a singlet quark distribution
and it is defined by

$\displaystyle
q_{\mathrm{S}}(x)=\sum_{i}q_{i}^{+}=\sum_{i}q_{i}(x)+\bar{q}_{i}(x),$ (22)

with $i$ is the quark flavor. The singlet quark distributions are expressed by

$\displaystyle\frac{\partial}{\partial\ln
Q^{2}}\begin{bmatrix}q_{\mathrm{S}}(x,Q^{2})\\\ g(x,Q^{2})\end{bmatrix}$
$\displaystyle=$
$\displaystyle\begin{bmatrix}\mathscr{P}_{qq}&\mathscr{P}_{qg}\\\
\mathscr{P}_{gq}&\mathscr{P}_{gg}\end{bmatrix}\otimes\begin{bmatrix}q_{\mathrm{S}}(x,Q^{2})\\\
g(x,Q^{2})\end{bmatrix}.$ (23)

Hence, the gluon distributions for the PS mesons can be obtained by solving
Eq. (23) numerically. In analogy to Taylor expansion series, the splitting
functions can be also expanded in terms of the running coupling constant,
$\alpha_{s}(Q^{2})$ in perturbative region and it then takes a form

$\displaystyle\mathscr{P}(z,Q^{2})$ $\displaystyle=$
$\displaystyle\left(\frac{\alpha}{2\pi}\right)\mathscr{P}^{(0)}(z)+\left(\frac{\alpha}{2\pi}\right)^{2}\mathscr{P}^{(1)}(z)+\cdot\cdot\cdot,$
(24)

where the first term in $\mathscr{P}^{(0)}(z)$ is the leading order (LO), the
second term in $\mathscr{P}^{(1)}(z)$ denotes the NLO. The NLO result is
expressed by

$\displaystyle\alpha_{s}$ $\displaystyle=$
$\displaystyle\frac{4\pi}{\beta_{0}}\frac{1}{\ln(Q_{\Lambda})}\left[1-\frac{\beta_{1}}{\beta_{0}}\frac{\ln\ln(Q_{\Lambda})}{\ln(Q_{\Lambda})}\right]+\mathscr{O}\left(\frac{1}{\ln^{2}(Q_{\Lambda})}\right),$
(25)

with

$\displaystyle Q_{\Lambda}$ $\displaystyle=$
$\displaystyle\frac{Q^{2}}{\Lambda_{QCD}^{2}},~{}~{}~{}~{}~{}\beta_{0}=\frac{11}{3}N_{c}-\frac{4}{3}N_{f},$
$\displaystyle\beta_{1}$ $\displaystyle=$
$\displaystyle\frac{34}{3}N_{c}^{2}-\frac{10}{3}N_{c}N_{f}-2C_{F}N_{f}.$ (26)

Here $N_{c}$, $N_{f}$ are the number of colors and the number of active
flavors, respectively and $C_{F}=\frac{4}{3}$. The value of $\Lambda_{QCD}$
depends on the number of active flavors and the renormalization scheme.

###  II.2 Nuclear matter NJL model

Here we present the equation of state (EoS) for SNM of the NJL model. Further
details can be found in Ref. [29] and references therein. Using the quark
bilinears that defined by $\bar{q}\Gamma
q=\langle\rho|\bar{q}\Gamma|\rho\rangle+:\bar{q}\Gamma q:$ with the vertices
$\Gamma=\mathds{1},\gamma^{\mu}$, the NJL effective Lagrangian in Eq. (1) is
modified in SNM as follows:

$\displaystyle\mathscr{L}_{\mathrm{SNM-NJL}}$ $\displaystyle=$
$\displaystyle\bar{q}(i\partial\\!\\!\\!/-M_{q}-V\\!\\!\\!\\!/)q-\frac{(M_{q}-m_{q})^{2}}{4G_{\pi}}+\frac{V_{\mu}V^{\mu}}{2G_{\omega}}$
(27) $\displaystyle+$ $\displaystyle\mathscr{L}_{I},$

where the in-medium constituent quark mass and the isoscalar-vector mean field
are defined respectively

$\displaystyle M_{q}$ $\displaystyle=$ $\displaystyle
m_{q}-2G_{\pi}\langle\rho|\bar{q}q|\rho\rangle,$ (28) $\displaystyle V^{\mu}$
$\displaystyle=$ $\displaystyle
2G_{\omega}\langle\rho|\bar{q}\gamma^{\mu}q|\rho\rangle=2G_{\omega}\langle
q^{\dagger}q\rangle,$ (29)

and $\mathscr{L}_{I}$ is the interaction Lagrangian. The $G_{\pi}$ and
$G_{\omega}$ are the scalar and vector coupling constants. The vector
potential is defined as $V=(V_{0},\mathbf{0})$.

Using the hadronization technique, the effective potential for SNM can be
obtained for the NJL Lagrangian. In the mean-field approach, the effective
potential is simply given by

$\displaystyle\mathscr{E}$ $\displaystyle=$
$\displaystyle\mathscr{E}_{V}-\frac{V_{0}^{2}}{4G_{\omega}}+4\int\frac{d^{3}p}{(2\pi)^{3}}\theta\left(p_{F}-|\mathbf{p}|\right)\epsilon_{p},$
(30)

where $\epsilon_{p}=\sqrt{M_{N}^{*2}+\mathbf{p}^{2}}+3V_{0}\equiv
E_{p}+3V_{0}$. Note that the effective nucleon mass as a function of the
constituent quark $M_{N}^{*}=M_{N}^{*}(M_{q}^{*})$ determined by solving the
quark-diquark bound state for the nucleon using the relativistic Faddeev
equation. The vacuum contribution of the quark loop in Eq. (30) is given by

$\displaystyle\mathscr{E}_{V}$ $\displaystyle=$ $\displaystyle
12i\int\frac{d^{4}k}{(2\pi)^{4}}\ln\left(\frac{k^{2}-M^{2}+i\epsilon}{k^{2}-M_{0}^{2}+i\epsilon}\right)+\frac{(M-m)^{2}}{4G_{\pi}}$
(31) $\displaystyle-$ $\displaystyle\frac{(M_{0}-m)^{2}}{4G_{\pi}}.$

In the present study, we consider the SNM at rest. Using the condition
$\partial\mathscr{E}/\partial V_{0}=0$, we then determine the value of $V_{0}$
as

$\displaystyle V_{0}$ $\displaystyle=$ $\displaystyle 6G_{\omega}\rho_{B},$
(32)

where $\rho_{B}=2p_{F}^{3}/3\pi^{2}$ is the baryon density with $p_{F}$ is the
Fermi momentum for the baryon. Also, the constituent quark mass $M_{q}$ for
fixed baryon density must satisfy the condition
$\frac{\partial\mathscr{E}}{\partial M_{q}}=0$ to give the similar expression
of the in-medium gap equation, as in Eq. (28).

A result for the effective nucleon and in-medium constituent quark masses as a
function of baryon density $\rho_{B}/\rho_{0}$ is shown in Fig. 2 (a). It
shows that the effective nucleon and in-medium quark constituent masses
decrease as the baryon density increases, as we expected, indicating the
partial restoration of the spontaneous breakdown of chiral symmetry in medium.

Thus the expression for the binding energy per nucleon is given by

$\displaystyle\frac{E_{B}}{A}$ $\displaystyle=$
$\displaystyle\frac{\mathscr{E}}{\rho_{B}}-M_{N0}.$ (33)

With the $M_{N0}$ is the nucleon mass in vacuum. A numerical result for
binding energy as a function of baryon density $\rho_{B}/\rho_{0}$ is depicted
in Fig. 2 (b). It shows that the SNM-NJL model can reproduce well the binding
energy per nucleon $E_{B}/A=$ -15.7 MeV at saturation density $\rho_{0}=$ 0.16
fm-3. The details of the obtained parameters in the SNM-NJL model will be
described in Sec. III.

r5cmt0.5cm(a) r5cmt0.5cm(b)

Figure 2:  (Color online) (a) The effective nucleon and constituent quark
masses as functions of the baryon density ($\rho_{B}/\rho_{0}$). (b) The
binding energy per nucleon as a function of the baryon density.

###  II.3 In-medium pion and kaon parton distributions

In this section the valence quark distributions of the PS mesons in SNM are
presented. Following the vacuum case, for general example, we present the
quark and antiquark distributions for the $K^{+}$-meson, which contains
[$u\bar{s}$] quarks, as

$\displaystyle u_{K^{+}}(x)$ $\displaystyle=$
$\displaystyle\frac{3\,g_{Kqq}^{*2}}{4\pi^{2}}\int d\tau\
e^{-\tau\left[x(x-1)\,m_{K}^{*2}+x\,M_{s}^{*2}+(1-x)\,M_{l}^{*2}\right]}$ (34)
$\displaystyle\mbox{}\times\left[\frac{1}{\tau}+x(1-x)\left[m_{K}^{*2}-(M_{l}^{*}-M_{s}^{*})^{2}\right]\right],$
$\displaystyle\bar{s}_{K^{+}}(x)$ $\displaystyle=$
$\displaystyle\frac{3\,g_{Kqq}^{*2}}{4\pi^{2}}\int d\tau\
e^{-\tau\left[x(x-1)\,m_{K}^{*2}+x\,M_{l}^{*2}+(1-x)\,M_{s}^{*2}\right]}$ (35)
$\displaystyle\mbox{}\times\left[\frac{1}{\tau}+x(1-x)\left[m_{K}^{*2}-(M_{l}^{*}-M_{s}^{*})^{2}\right]\right].$

The valence quark distributions of the $\pi^{+}$ in SNM can be easily obtained
by replacing $M_{s}^{*}\to M_{l}^{*}$ and $g_{Kqq}^{*}\to g_{\pi qq}^{*}$,
which gives the symmetry relation of $u_{\pi^{+}}(x)=\bar{d}_{\pi^{+}}(x)$.

The valence quark distributions of the PS mesons in SNM and vacuum with the
corresponding Bjorken variables $\tilde{x}_{a}$ and $x_{a}$, respectively, are
related by [20]

$\displaystyle q_{K^{+}}(x_{a})$
$\displaystyle=\frac{\epsilon_{F}}{E_{F}}q_{K^{+}}^{*}(\tilde{x}_{a})$ (36)

with

$\tilde{x}_{a}=\frac{\epsilon_{F}}{E_{F}}x_{a}-\frac{V^{0}}{E_{F}},$ (37)

where the quark energy in SNM
$\epsilon_{F}=\sqrt{k_{F}^{2}+M_{q}^{*2}}+V^{0}\equiv E_{F}+V^{0}$. The
$k_{F}$ is the quark Fermi momentum which is related to the quark density as
$\rho_{q}=2k_{F}^{3}/\pi^{2}$. For the antiquark, it is simply given by
$\epsilon_{F}=E_{F}-V^{0}$. The above formulas are valid only for light ($u$,
$d$) quarks in the present approach. Note that the in-medium valence quark
distributions satisfy the baryon number and momentum sum rules.

##  III Numerical result

Our numerical results for the in-medium and vacuum gluon and valence quark
distributions for the pion and kaon compared with the lattice QCD simulations
[15, 16] and the JAM phenomenology global fit QCD analysis [5] are presented.
The parameters of the present NJL model determined in vacuum are the coupling
constants $G_{\pi}$, $G_{\omega}$, $G_{\rho}$, $G_{a}$, $G_{s}$, and
$\Lambda_{\textrm{UV}}$ as in Refs. [27, 28]. In the PTR scheme, we fix the
regularization parameter of $\Lambda_{\textrm{IR}}=$ 0.240 GeV, which is in
order of the $\Lambda_{\textrm{QCD}}$, and opt the dressed light constituent
quark mass $M_{l}=$ 0.4 GeV. The remaining parameters are then fit to the
physical masses of the pion $m_{\pi}=$ 0.14 GeV, kaon $m_{K}=$ 0.495 GeV,
nucleon mass $M_{N}=M_{N0}=$ 0.94 GeV, and $\rho$-meson $m_{\rho}=$ 0.77 GeV
along with the decay constant of the pion $f_{\pi}=$ 0.093 GeV. This choice
gives $\Lambda_{\textrm{UV}}=$ 0.645 GeV, $G_{\pi}=$ 19.04 GeV-2, $G_{a}=$ 2.8
GeV-2, $G_{s}=$ 7.49 GeV-2, and $M_{qs}=$ 0.611 GeV, where $G_{s}$, $G_{a}$,
and $M_{\textrm{qs}}$ are respectively the scalar diquark coupling constant,
the axial-vector diquark coupling constant, and the strange constituent quark
mass. Note for the $\phi$ mass we get $m_{\phi}=$ 1.001 GeV, and the $G_{s}$
and $G_{a}$ coupling constants are determined in the Faddeev equations to fit
the vacuum nucleon mass $M_{N}$ and axial coupling constant of nucleon
$g_{A}=$ 1.267 [27]. Thus the scalar $M_{s}$ and axial $M_{a}$ diquark masses
are calculated using these model parameters and we obtain $M_{s}=$ 0.687 GeV
and $M_{a}=$ 1.027 GeV. The SNM parameters are determined to fit the binding
energy per nucleon $E_{B}/A=$ -15.7 MeV at saturation density $\rho_{0}=$ 0.16
fm-3, we then obtain $G_{\omega}=$ 6.03 GeV-2.

r5cmt0.5cm(a) r5cmt0.5cm(b) r5cmt0.5cm(a) r5cmt0.5cm(b)

Figure 3:  (Color online) (a) Valence quark and (b) gluon distributions of the
pion in SNM for various nuclear densities as well as in vacuum as functions of
the longitudinal momentum $x$ and the pion gluon distribution in vacuum
compared with the lattice QCD simulation [15] and JAM global fit analysis [5].
All DFs are evolved from the initial model scale of $Q_{0}^{2}=$ 0.16 GeV2 to
a scale $Q^{2}=$ 16 GeV2 using the NLO DGLAP QCD evolution. The experimental
data are the E615-Conway (blue shaded region) [31] and reanalysis E615 (orange
shaded region) [32]. Figure 4:  (Color online) The ratios of (a) the valence
quark distributions of the pion in SNM for various nuclear densities to the
pion vacuum valence-quark distributions and (b) for those of the pion gluon
distributions.

Results for the in-medium valence quark and gluon distributions of the pion as
a function of the light-front momentum fraction $x$ for various baryon
densities as well as that in vacuum that are evolved from the initial model
scale $Q_{0}^{2}=$ 0.16 GeV2 [28] to a scale at $Q^{2}=$ 16 GeV2 using the NLO
DGLAP evolution shown in Fig. 4 (a). It indicates that the valence quark
distribution for the pion seems less significant change in nuclear medium.
This result is consistent with the obtained result in Ref. [20]. The behaviors
of the in-medium valence quark PDFs are followed by the pion gluon
distributions in nuclear medium. In order to verify our approach for the gluon
PDFs, a comparison of our result for the pion gluon distribution in vacuum
with the lattice QCD simulation [15], and the JAM phenomenology global fit QCD
analysis [5] is presented in Fig. 4 (b). Our result shows a good agreement
with the results from the lattice QCD simulation (orange dotted line) [15],
and the JAM phenomenology global fit QCD analysis (blue dashed-dotted line)
[5].

r5cmt0.5cm(a) r5cmt0.5cm(b)
r5cmt0.5cm(c)

Figure 5:  (Color online) (a) Up valence-quark distributions of the kaon, (b)
the strange valence-quark distributions of the kaon, and (c) the gluon
distributions of the kaon.

The ratios of the pion valence-quark distributions in SNM to those in vacuum
for various densities are shown in Fig. 4 (a). Figure 4 (a) indicates that the
ratios decrease as the nuclear density increases up to $x\simeq 0.2$, which is
consistent with the finding in Ref. [30]. However, the reduction is relatively
small for $x\lesssim 0.2$. This is because, in the region of $x\lesssim 0.2$,
the pion gluon distribution is expected relatively large and it becomes
growing in nuclear medium that might be suppressed the valence quark DFs in
the pion. In contrast, for $x\gtrsim$ 0.2, the ratios increase with respect to
the nuclear matter density. The ratios of the gluon distributions for the pion
are depicted in Fig. 4 (b) in the same manner. It shows that the ratios of the
pion gluon distributions in nuclear medium for various densities to those in
vacuum increase relatively weakly with respect to the density.

For the kaon case, the kaon valence-quark distributions are evolved from the
initial model scale $Q_{0}^{2}=$ 0.16 GeV2 to $Q^{2}=$ 16 GeV2 using the NLO
DGLAP evolution and depicted in Fig. 5 (a). We find that the up valence-quark
distributions for the kaon increase with respect to the density. The
enhancement of the up valence-quark distributions for the kaon in nuclear
medium are relatively significant. It is more pronounced at higher densities.
The results for the strange quark distributions in nuclear medium for various
densities as well as those in vacuum after evolved at $Q^{2}=$ 16 GeV2 are
shown in Fig. 5 (b). We find that the strange valence-quark distributions of
the kaon decrease as the density increases in the range $0\lesssim x\lesssim
0.4$. Thus, at $x\gtrsim 0.4$, the strange valence-quark distributions of the
kaon begin to increase as the density increases. The gluon distributions of
the kaon in nuclear medium are shown in Fig. 5 (c). In contrast to the gluon
contents of the pion in nuclear medium for various densities, it seems that
the gluon distributions of the kaon in nuclear medium increase more
significantly as the density increases. The enhancement on the gluon
distributions for the kaon in nuclear medium can be clearly seen in Fig. 6
(a). Also, the ratios of the up valence-quark distributions of the kaon to
those of the pion in nuclear medium are given in Fig. 6 (b). It shows that the
ratios of the up valence-quark distributions of the kaon to those of the pion
in nuclear medium increase as the density increases and it decreases as the
$x$ increases.

r5cmt0.5cm(a) r5cmt0.5cm(b)
r5cmt0.5cm(c) r5cmt0.5cm(d)

Figure 6:  (Color online) The ratios for (a) the gluon distributions of the
kaon in SNM to the vacuum gluon distributions, (b) the up valence-quark
distributions of the kaon to the up valence-quark distributions of the pion,
(c) the gluon distributions of the kaon to the gluon distributions of the
pions and (d) the gluon distribution differences for the kaon and pions in SNM
for various nuclear densities.

The ratios of the gluon distributions of the kaon to those of the pion as a
function of $x$ for various nuclear densities are shown in Fig. 6 (c). It can
be seen that the ratios of the gluon distributions in the kaon increase with
respect to the densities. However, it decreases as the $x$ increases. It
implies the gluon distributions in the pion are relatively larger than those
in the kaon. This can be understood because the strange quark, which has
heavier mass than the light quark, radiates less instantly than the light
quark and it radiates soft gluons than the light quark does. This finding is
compatible with the result obtained in the Dyson-Schwinger equations (DSEs)
model [8], which also pointed out in Ref. [10]. Figure 6 (c) depicts that the
larger gluon distribution contribution in the pion occurs at
$\rho_{B}\lesssim\rho_{0}$, which gives the ratios of
$g_{K}(x,\rho_{B})/g_{\pi}(x,\rho_{B})\lesssim$ 1\. Unexpectedly, we find that
the ratios of $g_{K}(x,\rho_{B})/g_{\pi}(x,\rho_{B})\gtrsim$ 1 for
$\rho_{B}>\rho_{0}$. This indicates that the gluon distribution in the kaon
becomes larger than those in the pion for the higher density region
($\rho_{B}>\rho_{0}$).

To more clearly understand this unexpected gluon contribution, we calculate
the gluon distribution differences for the pion and kaon as shown in Fig. 6
(d). We find that Figure 6 (d) supports what we found in Fig. 6 (c). At
$\rho_{B}\lesssim\rho_{0}$, it gives $g_{K}(x,\rho_{B})-g_{\pi}(x,\rho_{B})<$
0 (negative value), but for $\rho_{B}>\rho_{0}$, we find that
$g_{K}(x,\rho_{B})-g_{\pi}(x,\rho_{B})>$ 0 (positive value). The negative
values of the $g_{K}(x,\rho_{B})-g_{\pi}(x,\rho_{B})$ clearly explains that
the gluon content in the pion is significantly larger than that in the kaon in
vacuum and the nuclear density up to the saturation density,
$\rho_{B}\simeq\rho_{0}$. The positive values mean that the gluon content in
the pion is relatively lower than that in the kaon in vacuum and the nuclear
medium at $\rho_{B}>\rho_{0}$.

##  IV Summary

To summarize, we have studied the in-medium and vacuum valence-quark and gluon
distributions of the pion and kaon in the Nambu–Jona-Lasinio model with the
help of the proper-time regularization scheme, simulating a QCD confinement.
To describe nuclear matter, we adopt the SNM-NJL model. We then computed the
valence quark distributions of the pion and kaon in nuclear medium for various
nuclear densities as well as in vacuum. In this work, at initial scale
$Q_{0}^{2}$, the gluon distributions in the pion and kaon are taken to be
zero, since there is no gluon dynamics in the NJL model. So, the gluon
distributions in the pion and kaon in SNM for various nuclear densities as
well as in vacuum are purely generated from the NLO DGLAP QCD evolution.

We have determined the gluon and quark distributions in SNM at a higher scale
by evolving them from the initial scale $Q_{0}^{2}=$ 0.16 GeV2 to a
factorization scale at $Q^{2}=$ 16 GeV2. Our result for the vacuum gluon and
valence-quark distributions of the pion have good agreement with the results
obtained from the JAM phenomenology global fit QCD analysis [5], the lattice
QCD simulation [15], and the experimental data [31], respectively. The results
for the in-medium gluon and valence-quark distributions in the pion for
various nuclear densities shows that their distributions are almost unchanged
in nuclear medium.

The results for the kaon case, the gluon and up valence-quark distributions
for various nuclear densities significantly modify in nuclear medium, in
particular, for the in-medium up valence-quark distributions in the kaon. The
in-medium strange valence-quark distributions in the kaon for various
densities decrease as the nuclear densities increase at $x\lesssim 0.4$. In
contrast, it becomes increasing as the nuclear density increases at $x\gtrsim
0.4$.

Our results clearly shows that the vacuum gluon distributions in the kaon is
smaller than those in the pion, which is consistent with the DSEs result in
Ref. [8]. Also, we found that the in-medium gluon distributions in the kaon is
smaller than those in the pion up to $\rho_{B}\simeq\rho_{0}$. Surprisingly,
for higher nuclear densities $\rho_{B}>\rho_{0}$, the gluon distributions in
the kaon are found to be bigger than those in the pion as shown in Fig. 6 (d).
This finding may be useful and relevant information for the compact star
(neutron star) or heavy-ion collisions studies. In such cases, the NJL model
serves as a very useful tool to guide possible future computations for the
lattice QCD and other more sophisticated theoretical approaches.

###### Acknowledgements.

P. T. P. H thanks X. G. Wang and Ian Cloët for useful discussions. This work
was supported by the National Research Foundation of Korea (NRF) funded by the
Korea government (MSIT) Grants No. 2018R1A5A1025563 and No. 2019R1A2C1005697.

## References

  * [1] E. L. Berger and S. J. Brodsky,  Quark Structure Functions of Mesons and the Drell-Yan Process,  Phys. Rev. Lett. 42, 940-944 (1979).
  * [2] G. Xie, C. Han, R. Wang and X. Chen,  Tackling the kaon structure function at EicC,  [arXiv:2109.08483 [hep-ph]].
  * [3] J. M. M. Chavez, V. Bertone, F. D. Borrero, M. Defurne, C. Mezrag, H. Moutarde, J. Rodríguez-Quintero and J. Segovia,  Pion GPDs: A path toward phenomenology,  [arXiv:2110.06052 [hep-ph]].
  * [4] C. S. R. Costa, A. Freese, I. C. Cloët, B. El-Bennich, G. Krein and P. C. Tandy,  Intrinsic glue and Wilson lines within dressed quarks,  Phys. Rev. C 104, no.4, 045201 (2021).
  * [5] P. C. Barry, N. Sato, W. Melnitchouk and C. R. Ji,  First Monte Carlo Global QCD Analysis of Pion Parton Distributions,  Phys. Rev. Lett. 121, no.15, 152001 (2018).
  * [6] J. Arrington, C. A. Gayoso, P. C. Barry, V. Berdnikov, D. Binosi, L. Chang, M. Diefenthaler, M. Ding, R. Ent and T. Frederico, et al. Revealing the structure of light pseudoscalar mesons at the electron–ion collider,  J. Phys. G 48, no.7, 075106 (2021).
  * [7] S. i. Nam,  Parton-distribution functions for the pion and kaon in the gauge-invariant nonlocal chiral-quark model,  Phys. Rev. D 86, 074005 (2012).
  * [8] C. Chen, L. Chang, C. D. Roberts, S. Wan and H. S. Zong,  Valence-quark distribution functions in the kaon and pion,  Phys. Rev. D 93, no.7, 074021 (2016).
  * [9] P. T. P. Hutauruk, W. Bentz, I. C. Cloët and A. W. Thomas,  Charge Symmetry Breaking Effects in Pion and Kaon Structure,  Phys. Rev. C 97, no.5, 055210 (2018).
  * [10] M. Ding, K. Raya, D. Binosi, L. Chang, C. D. Roberts and S. M. Schmidt,  Symmetry, symmetry breaking, and pion parton distributions,  Phys. Rev. D 101, no.5, 054014 (2020).
  * [11] A. C. Aguilar, Z. Ahmed, C. Aidala, S. Ali, V. Andrieux, J. Arrington, A. Bashir, V. Berdnikov, D. Binosi and L. Chang, et al. Pion and Kaon Structure at the Electron-Ion Collider,  Eur. Phys. J. A 55, no.10, 190 (2019).
  * [12] D. P. Anderle, V. Bertone, X. Cao, L. Chang, N. Chang, G. Chen, X. Chen, Z. Chen, Z. Cui and L. Dai, et al. Electron-ion collider in China,  Front. Phys. (Beijing) 16, no.6, 64701 (2021).
  * [13] B. Adams, C. A. Aidala, R. Akhunzyanov, G. D. Alexeev, M. G. Alexeev, A. Amoroso, V. Andrieux, N. V. Anfimov, V. Anosov and A. Antoshkin, et al. Letter of Intent: A New QCD facility at the M2 beam line of the CERN SPS (COMPASS++/AMBER),  [arXiv:1808.00848 [hep-ex]].
  * [14] J. Lan, K. Fu, C. Mondal, X. Zhao and j. P. Vary,  Light mesons with one dynamical gluon on the light front,  [arXiv:2106.04954 [hep-ph]].
  * [15] Z. Fan and H. W. Lin,  Gluon Parton Distribution of the Pion and Nucleon from Lattice QCD,  [arXiv:2110.14471 [hep-lat]].
  * [16] I. Novikov, H. Abdolmaleki, D. Britzger, A. Cooper-Sarkar, F. Giuli, A. Glazov, A. Kusina, A. Luszczak, F. Olness and P. Starovoitov, et al. Parton Distribution Functions of the Charged Pion Within The xFitter Framework,  Phys. Rev. D 102, no.1, 014040 (2020).
  * [17] R. J. Holt and C. D. Roberts,  Distribution Functions of the Nucleon and Pion in the Valence Region,  Rev. Mod. Phys. 82, 2991-3044 (2010).
  * [18] R. D. Ball, E. R. Nocera and J. Rojo,  The asymptotic behaviour of parton distributions at small and large $x$,  Eur. Phys. J. C 76, no.7, 383 (2016).
  * [19] J. P. B. C. de Melo, K. Tsushima, B. El-Bennich, E. Rojas and T. Frederico,  Pion structure in the nuclear medium,  Phys. Rev. C 90, no.3, 035201 (2014).
  * [20] P. T. P. Hutauruk, J. J. Cobos-Martínez, Y. Oh and K. Tsushima,  Valence-quark distributions of pions and kaons in a nuclear medium,  Phys. Rev. D 100, no.9, 094011 (2019).
  * [21] X. G. Wang, W. Bentz, I. C. Cloët and A. W. Thomas,  Gluon EMC Effects in Nuclear Matter,  [arXiv:2109.03591 [hep-ph]].
  * [22] M. Miyama and S. Kumano,  Numerical solution of Q**2 evolution equations in a brute force method,  Comput. Phys. Commun. 94, 185-215 (1996).
  * [23] Y. Ninomiya, W. Bentz and I. C. Cloët,  Transverse-momentum-dependent quark distribution functions of spin-one targets: Formalism and covariant calculations,  Phys. Rev. C 96, no.4, 045206 (2017).
  * [24] H. H. Matevosyan, A. W. Thomas and W. Bentz,  Dihadron fragmentation functions within the Nambu–Jona-Lasinio–jet model,  Phys. Rev. D 88, no.9, 094022 (2013).
  * [25] P. T. P. Hutauruk, Y. Oh and K. Tsushima,  Electroweak properties of pions in a nuclear medium,  Phys. Rev. C 99, no.1, 015202 (2019).
  * [26] P. T. P. Hutauruk and K. Tsushima,  Electroweak properties of kaons in a nuclear medium,  [arXiv:1910.08133 [nucl-th]].
  * [27] T. Tanimoto, W. Bentz and I. C. Cloët,  Massive Neutron Stars with a Color Superconducting Quark Matter Core,  Phys. Rev. C 101, no.5, 055204 (2020).
  * [28] P. T. P. Hutauruk, I. C. Cloet and A. W. Thomas,  Flavor dependence of the pion and kaon form factors and parton distribution functions,  Phys. Rev. C 94, no.3, 035201 (2016).
  * [29] W. Bentz and A. W. Thomas,  The Stability of nuclear matter in the Nambu-Jona-Lasinio model,  Nucl. Phys. A 696, 138-172 (2001).
  * [30] L. Chang and C. D. Roberts,  Regarding the distribution of glue in the pion,  Chin. Phys. Lett. 38, no.8, 081101 (2021).
  * [31] J. S. Conway, C. E. Adolphsen, J. P. Alexander, K. J. Anderson, J. G. Heinrich, J. E. Pilcher, A. Possoz, E. I. Rosenberg, C. Biino and J. F. Greenhalgh, et al. Experimental Study of Muon Pairs Produced by 252-GeV Pions on Tungsten,  Phys. Rev. D 39, 92-122 (1989).
  * [32] M. Aicher, A. Schafer and W. Vogelsang,  Soft-gluon resummation and the valence parton distribution function of the pion,  Phys. Rev. Lett. 105, 252003 (2010).