# HEROS: Hierarchical Exploration with Online Subregion Updating for 3D
Environment Coverage

Shijun Long, Ying Li, Chenming Wu, Bin Xu, and Wei Fan This work was supported
in part by the National Natural Science Foundation of China under Grant
52102449, in part by the Natural Science Foundation of Chongqing under Grant
CSTB2023NSCQ-MSX0550, and in part by the Beijing Institute of Technology
Research Fund Program for Young Scholars and S&T Program of Hebei under Grant
21567606H.S. Long, Y. Li, B. Xu, and F. Wei are with the School of Mechanical
Engineering, Beijing Institute of Technology, Beijing, China. {sj_long,
ying.li, bitxubin<EMAIL_ADDRESS>Wu is with RAL, Baidu Research.
<EMAIL_ADDRESS>
###### Abstract

We present an autonomous exploration system for efficient coverage of unknown
environments. First, a rapid environment preprocessing method is introduced to
provide environmental information for subsequent exploration planning. Then,
the whole exploration space is divided into multiple subregion cells, each
with varying levels of detail. The subregion cells are capable of
decomposition and updating online, effectively characterizing dynamic unknown
regions with variable resolution. Finally, the hierarchical planning strategy
treats subregions as basic planning units and computes an efficient global
coverage path. Guided by the global path, the local path that sequentially
visits the viewpoint set is refined to provide an executable path for the
robot. This hierarchical planning from coarse to fine steps reduces the
complexity of the planning scheme while improving exploration efficiency. The
proposed method is compared with state-of-art methods in benchmark
environments. Our approach demonstrates superior efficiency in completing
exploration while using lower computational resources.

###### Index Terms:

Search and rescue robots, mapping, autonomous agents.

##  I Introduction

Autonomous exploration using ground mobile robots has received considerable
attention in recent years for its wide applications, including searching and
rescuing, environmental monitoring, and industrial automation. The ability of
robots to navigate and explore unknown environments without human intervention
enhances operational efficiency and ensures personnel safety in hazardous
conditions. Although several autonomous exploration methods have been proposed
over the decades [1], [2], [3], the task remains challenging in large-scale
and topologically complex scenarios. In such scenarios, the problem becomes
computationally complex, and the robot needs to navigate between areas in a
limited time while avoiding redundant exploration.

The mainstream exploration frameworks mainly include frontier-based and
sampling-based methods. In autonomous exploration, the frontier is defined as
the boundary between free and unknown areas [1]. These two approaches guide
the robot to cover the environment by continuously searching for frontiers in
the map or randomly sampling viewpoints and evaluating their potential
revenue. However, several imitations still hinder the development of
autonomous exploration:

High Computational Cost. Both frontier detection and revenue evaluation are
computationally expensive, requiring frequent checks of each voxel’s occupancy
status or performing ray tracing [4, 2, 5]. This presents a substantial
challenge for ground-based platforms equipped with edge computing devices to
meet real-time requirements.

Low Efficiency. Most methods [2, 6, 1] employ a greedy strategy of direct
traveling to the frontier or viewpoint with the highest gain, disregarding the
global coverage path. However, it is crucial to consider the potential impact
on future exploration planning when selecting the next target, to avoid
inefficiencies caused by shortsighted decisions. While some approaches [7],
[8] employ the Travelling Salesman Problem (TSP) to obtain the global coverage
path, as the number of viewpoints increases, the TSP becomes larger and
heavier to solve due to its NP-Complete natural [5].

Figure 1: The real-world experiment conducted in an underground garage using
real ground platforms. The yellow dot denotes the start position, and the
white line indicates the robot’s trajectory.

Taking inspiration from the aforementioned concerns, we propose an efficient
hierarchical exploration planner to enable robots to explore complex unknown
environments autonomously. First, we design a rapid environmental
preprocessing method that incrementally extracts frontiers, continuously
constructs the sparse graph of the environment, and samples high-quality
viewpoints. To effectively represent the continuously changing unknown areas
during the exploration process, we adopt a subregion-based regional division
approach. This approach decomposes the entire unknown environment into several
subregion cells with a hierarchical structure, which can be further subdivided
to adapt to the changing unknown areas. Additionally, the Subregion
Information Structure (SIS) is introduced to record basic information within
the subregion space, facilitating subsequent exploration planning. Finally,
the subregion is adopted as the basic planning unit, and the hierarchical
planning strategy is used for efficient exploration. The hierarchical planning
includes global coverage planning and local planning. The global coverage path
is the shortest travel path that enters each subregion in sequence, while the
local path is the viewpoint access path refined under the guidance of the
global path. The above process is continuously executed until the entire
environment is fully explored.

The main contributions of this paper are as follows:

  * •

A subregion-based regional division method with a hierarchical structure is
proposed, allowing for adaptive decomposition of the subregion to capture the
dynamic changes of unknown regions.

  * •

A hierarchical planning strategy is adopted, achieving a significant balance
between coverage exploration efficiency and computational complexity.

  * •

Extensive simulation and real-world experiments are conducted to verify the
theoretical and practical feasibility of the proposed method. Our approach
demonstrates the fastest exploration with minimal computational resource
consumption across diverse scenarios.

The source code will be released for further research111Available at
https://github.com/bit-lsj/HEROS.git.

##  II Related Works

The field of autonomous robotic exploration has seen significant advancements
in recent years, with a variety of approaches being proposed to tackle the
challenges of unknown and complex environments. Early work in this domain
focused on heuristic-based methods, such as frontier-based exploration [1, 3,
9]. This type of method typically selects the nearest frontier or the frontier
with the highest information gain as the next target. They have worked in
simple environments but often fall short in complex large-scale scenes due to
their greediness. Besides, these methods still require a significant amount of
frontier detection and clustering, making it difficult to meet real-time
requirements for unmanned platforms. [10] introduces an incremental frontier
detection method and performs global coverage path planning after clustering
frontiers. Nevertheless, in large-scale environments, the algorithm’s runtime
speed is severely slowed down due to a large number of clusters and solving
TSP. [11] greatly improves planning efficiency by focusing on a specific range
of viewpoints for exploration planning. However, it lacks global path guidance
and can still lead to inefficient exploration in some cases.

The sampling-based method evaluates sampled viewpoints using a utility
function and selects the viewpoint or branch with the highest gain to observe
the unknown region. The idea of the next best views [12] is first introduced
into autonomous exploration by [2]. This method utilizes rapidly random trees
(RRT) to explore 3D spaces, selecting the branch with the highest gain for
execution. This method has been continuously improved by [4], [5], but they
all suffer from the trouble of expensive computational costs. The majority of
computational resources for these methods are consumed in evaluating the gains
of viewpoints or branches, which requires estimating the quantity of a large
number of unknown voxels.

Several methods [13], [14], [15], [16] are proposed for the concepts of
regional division and hierarchical planning. However, these methods are either
designed for task allocation in multi-robot collaborative exploration or
utilize fixed-sized regions that cannot adapt to dynamic changes in unknown
areas. In [15], the global path is introduced to provide a more comprehensive
perspective for exploration. However, the greedy strategy is still adopted in
the local planning stage, which can still lead to short-sighted behavior and
thus reduce exploration efficiency. The concept of region hierarchical
division in this paper is inspired by [14], which utilizes hgrid space for
muti-UAV interaction and task assignment. Our approach of the regional
division is aimed at serving the global coverage path planning, which is part
of the hierarchical planning.

Learning-based techniques have also emerged as promising solutions for
autonomous exploration [17], [18], [19]. These methods leverage data-driven
models to predict optimal exploration strategies but often lack robustness and
adaptability to unfamiliar scenarios. This leads to performance deterioration
in complex and unfamiliar environments.

Figure 2: System overview. The environmental preprocessing is first conducted,
which includes frontier detection, building sparse graph, and viewpoint
generation. Then, the entire exploration environment is divided into several
subregion cells, and the SIS of each cell is stored for global coverage path
planning and local path optimization. Finally, the robot follows the local
path to explore.

##  III System Overview.

The overview of the proposed exploration framework is presented in Fig. 2. The
exploration task is decomposed into several sub-tasks, which are then
organized into a timeline for sequential execution. The task starts with rapid
environment preprocessing (Sec. IV), which includes frontier extraction,
sparse graph construction, and viewpoint generation. The unknown space is then
divided into hierarchical subregion cells with recorded environmental
information (Sec. V). Using these cells, a global coverage path is generated,
guiding the refinement of the local path to cover each viewpoint efficiently
(Sec. VI). The exploration ends when no viewpoints remain.

##  IV Rapid Environmental Preprocessing

This section introduces a real-time and scalable preprocessing method for
large-scale scenarios to support efficient exploration planning. It covers
three main components: environment modeling and frontier detection, sparse
graph construction, and viewpoint generation.

###  IV-A Mapping and Frontier Detection

In this work, the UFOMap [20] is utilized for environmental representation,
for its advantages of low memory consumption and fast traversal speed.
Additionally, with the use of UFOMap, it becomes convenient to fetch the
boundary of the updated region in the map and perform collision checking using
rapid ray tracing.

As in state-of-art works [10, 21], an incremental frontier detection method is
employed in this work to avoid processing the entire map. Specifically, a set
$\mathcal{C}_{v}$ of all voxels for which the occupation value has changed is
obtained from the UFOMap, and the state of each voxel $V\in\mathcal{C}_{v}$ is
checked. If a voxel $V$ is in the free state (denoted as $V_{free}$) and an
adjacent voxel is in the unknown state (denoted as $V_{un}$) at the same
height, then $V$ can be considered as a frontier voxel. During the exploration
process, each frontier voxel located in the map updated region
$\mathcal{B}_{u}$ is checked to determine if it can still be a frontier, and
those invalid frontiers are removed. This method only requires checking a
small number of voxels in the map, significantly reducing the workload for
frontier extraction.

###  IV-B Dynamic Sparse Graph

Given that our exploration method requires extensive path search operations,
which are time-consuming in the 3D voxel map. Thus, a dynamic sparse graph
$\mathcal{G}=(N,E)$ is maintained as a sparse representation of the
environment. Here, $N$ represents the nodes in the graph, and $E$ denotes the
undirected weighted edge connecting adjacent nodes. During the exploration,
candidate nodes $N_{c}$ are generated through the random biased sampling
method based on the terrain map $\mathcal{A}$. The sampled node that satisfies
the conditions in Eq. 1 will undergo collision checking, and then be added to
$\mathcal{G}$. The weight of the edge $E$ is set as the Euclidean distance
between two nodes. All the nodes are stored in R*-tree [22] to speed up the
computational efficiency. It is also important to perform collision checks for
all edges within the boundary $\mathcal{B}_{u}$ during the exploration, and
remove those that become invalid due to dynamic changes in the environment.

$\begin{cases}\textit{$N_{c}$ is in Accessible}(\mathcal{A})\\\
\gamma<D_{E}(N_{\text{nearest}},N_{c})<\delta\\\ \end{cases}$ (1)

where ${Accessible}(\mathcal{A})$ denotes the accessible region in
$\mathcal{A}$, since all nodes should be reachable by the robot. The
parameters $\delta$ and $\gamma$ represent the minimum and maximum allowable
distances between $N_{c}$ and the nearest node $N_{\text{nearest}}$ in
$\mathcal{G}$, respectively. This can prevent the graph from becoming over
dense or sparse.

Figure 3: Visualization effect of environmental preprocessing. During the
exploration process, new nodes are continuously sampled to expand the sparse
graph $\mathcal{G}$. Viewpoints are obtained by sampling from $\mathcal{G}$
through generating collision-free spheres. All frontier points within the
sphere can be observed from the viewpoint.

###  IV-C Viewpoint Generation

In most approaches [2, 5, 10], the viewpoint is utilized to observe unknown
voxels within the visible range of sensors, thus continuously extending the
frontier. However, this process requires frequent and expensive ray tracing to
assess the visibility of the unknown voxels from the viewpoint. Inspired by
[23, 24], we employ the collision-free sphere to generate high-quality
viewpoints to observe the unknown environment. This approach aims to reduce
the number of viewpoints, thus simplifying the complexity of subsequent
exploration planning. The center $\mathbf{p}_{s}$ of the collision-free sphere
is located at the frontier voxel, with its radius defined as:

$r=\min(\|\mathbf{p}_{s}-\mathbf{p}_{o}\|_{2},\frac{r_{\text{sensor}}}{2})$
(2)

where $\mathbf{p}_{o}$ is the position of the nearest occupied voxel, and
$r_{\text{sensor}}$ is the maximum scanning radius of the LiDAR sensor. Due to
the absence of obstacles within the sphere and its convex nature, any frontier
voxel in the sphere can be observed from a viewpoint inside. Thus, performing
expensive ray tracing on every frontier voxel to determine its visibility from
the viewpoint is unnecessary. The overall generation procedure of the
collision-free sphere is similar to [24]. Then, the node $N$ located inside
the sphere and closest to its center is selected as the viewpoint. This search
process is achieved by K-Nearest Neighbor search in $\mathcal{G}$. Fig. 3
shows the dynamic sparse graph and the viewpoints generated based on
collision-free spheres.

##  V Subregion-based Regional Division

Instead of directly solving the visit order of viewpoints or frontiers, we use
subregion cells as basic planning units for global coverage. Each candidate
cell, with varying levels of detail, is managed through online decomposition
and updating for a multi-resolution representation of the unknown space.

Figure 4: Schematic representation of the hierarchy of subregion cells.
Subregion cells can be further subdivided to represent unknown regions using
smaller resolutions up to the maximum hierarchy $l_{\text{max}}$. Figure 5: An
example of the regional division with two levels. The original partitioned
regions will be further subdivided if the known voxel ratio inside exceeds a
threshold. Regions with no viewpoints inside will be considered as being in
inactive status.

###  V-A Subregion hierarchical structure

We utilize a quadtree-like form to continuously subdivide the whole
environment into multiple subregions, characterizing the unknown regions
through different levels of detail. Similar to the octree’s structural
division, a subregion cell at level $l$ can be recursively subdivided into 2 ×
2 subcells in level $l+1$. This method aids in representing the area to be
explored with varying resolutions. Fig. 4 shows the hierarchical structure of
the subregion cell.

At the beginning of the exploration, the total exploration space $\mathcal{O}$
is discretized into $S_{o}=S_{o,x}\times S_{o,y}$ cells along the axis-aligned
bounding box (AABB) of the map. For any candidate subregion cell
$S_{i}\in\left\\{S_{o}\right\\}$, as the exploration process proceeds, the
volume of the unknown regions in $S_{i}$ gradually decreases. As a
consequence, the original divided cell fails to accurately represent the
unknown regions. Thus, the $S_{i}$ is further subdivided into four subcells to
describe the unknown regions at a smaller resolution. The recursive
subdivision method for subregions can describe unknown areas at low cost and
decompose online as the exploration process progresses. This provides
sufficient flexibility for representing unknown regions. Fig. 5 depicts an
example of the regional division.

Note that the subregion cells are only rough representations of unknown
regions, which are used as basic exploration units to achieve global coverage.
If the minimum area of a subregion cell is too small, a large number of
subregion cells will be formed during the exploration process. This increases
the complexity of the subsequent global path solving. Let $W_{\text{min}}$ and
$W_{\text{max}}$ denote the minimum and maximum side lengths of the subregion
cell, respectively. During the exploration, $W_{\text{min}}$ and
$W_{\text{max}}$ should be set in an appropriate range, and the maximum level
$l_{\text{max}}$ can be calculated as:

$l_{\text{max}}=\left\lceil\log_{2}\left(\frac{W_{\text{max}}}{W_{\text{min}}}\right)+1\right\rceil$
(3)

where $\lceil\cdot\rceil$ denotes the ceiling function.

###  V-B Subregion Information Structure

The Subregion Information Structure $\mathbf{SI_{i}}$ is composed of elemental
information contained in each subregion cell $S_{i}$. The hierarchy level
$l_{i}$ and axis-aligned bounding box $\mathbf{B}_{i}$ are recorded when
initializing $S_{i}$. The proportion of known voxels
$\mathbf{R}_{\text{known},i}$ and the viewpoint average position
$\mathbf{C}_{\text{avg},i}$ of $S_{i}$ are also computed to locate the area to
be explored. To serve subsequent global coverage path planning, the traveling
cost $\mathbf{L}_{\text{cost},i}$ between $S_{i}$ and other subregions is also
computed. The viewpoints $\mathbf{VP}_{i}$ within $S_{i}$ are stored in
$\mathbf{SI_{i}}$ for subsequent local path optimization. In addition, each
subregion cell holds a status $stat_{i}$ between “active” and “inactive”.
Subregion cells containing at least one viewpoint are categorized as active,
while those without viewpoints are considered inactive. Only subregions with
active status are eligible for inclusion in the global coverage path planning
process. Table. I presents the meaning of each data information in SIS.

TABLE I: Information contained by SIS $\mathbf{SI_{i}}$ of subregion cell
$S_{i}$.

Information | Note
---|---
$l_{i}$ | Hierarchical level of the cell
$\mathbf{B}_{i}$ | Axis-aligned bounding box of $S_{i}$
$\mathbf{R}_{\text{known},i}$ | Proportion of known voxels in the subregion
$\mathbf{C}_{\text{avg},i}$ | Average position of viewpoints within $S_{i}$
$\mathbf{L}_{\text{cost},i}$ | Traveling cost from $\mathbf{C}_{\text{avg},i}$ of $S_{i}$ to other cells
$\mathbf{VP}_{i}$ | Viewpoints in subregion
$stat_{i}$ | ”active” or ”inactive” status of $S_{i}$

###  V-C Online Subregion Updating

We maintain a set $\mathcal{V}$ that stores all divided subregions. At the
start of the exploration, the entire unknown space $\mathcal{O}$ is divided
into $S_{o}=S_{o,x}\times S_{o,y}$ cells at level $l=1$. During the
exploration process, as the map is updated, the $\mathbf{SI_{i}}$ of $S_{i}$
also changes accordingly. Therefore, all cells intersecting with the map
update region $\mathcal{B}_{u}$ need to recompute their SIS. Algorithm 1
outlines the online subregion cell updating process. First, traverse the set
$\mathcal{V}$ to filter the cell $S_{i}$ that overlaps with $\mathcal{B}_{u}$
(Line 2). Then, some information such as the proportion of known voxels
$\mathbf{R}_{\text{known},i}$ and the status $stat_{i}$ is updated (Line 4).
For $S_{i}$ in active status, if $\mathbf{R}_{\text{known},i}$ exceeds a
certain threshold $\sigma_{k}$ and is not at $l_{\text{max}}$, it is
subdivided into four subcells (Line 8). Then, the $S_{i}$ is removed from
$\mathcal{V}$, and the SIS of the subcells is calculated (Lines 9-12).
Otherwise, the passable path from $S_{i}$ to other subregion cells is searched
in $\mathcal{G}$ to update traveling cost $\mathbf{L}_{\text{cost},i}$ (Line
14). Note that the update of $\mathbf{L}_{\text{cost}}$ between subregion pair
is mutual, so the traveling cost from another subregion $S_{j}$ to $S_{i}$
also needs to be updated, and
$\mathbf{L}_{\text{cost},i,j}=\mathbf{L}_{\text{cost},j,i}$.

For $N_{\text{act}}$ subregions in active status, if we calculate the
traveling cost $\mathbf{L}_{\text{cost}}$ between all subregion pairs,
$O(N_{\text{act}}^{2})$ A* path searches are required at each iteration. In
most cases, the number of $N_{\text{act}}$ is not large. However, in extensive
environments, this can still be very computationally intensive and may cause
the robot to freeze if the subregion pairs are far apart. Theoretically, only
the subregion cells intersecting the updated region $\mathcal{B}_{u}$ are
affected, so only the cost data of these cells need to be updated. Assume that
$M$ subregion cells are affected, in this case, only $O(MN_{\text{act}})$ path
searches are performed. Practically, the number of $M$ is small, resulting in
a small number of path search computations.

// Initialize $\mathcal{V}=\left\\{S_{o}\right\\}$ at the beginning

1 for  _$S_{i}\in\mathcal{V}$_ do

2 if  _$\neg$ isOverlap($S_{i},\mathcal{B}_{u}$)_ then

3 continue

4 end if

5 updateInfo($S_{i}$)

6 if  _$\neg$ isActive($S_{i}$)_ then

7 continue

8 end if

9 if  _$\mathbf{R}_{\text{known},i} >\sigma_{k}\land l_{i}<l_{\text{max}}$_
then

10 $S_{\text{div}}\leftarrow$ subdivideCell($S_{i}$)

11 $\mathcal{V}$.erase($S_{i}$)

12 for  _$S_{l}\in S_{\text{div}}$_ do

13 calculateSIS($S_{l}$)

14 $\mathcal{V}$.append($S_{l}$)

15

16 end for

17

18 end if

19 else

20 updateDist($S_{i}$)

21

22 end if

23

24 end for

Algorithm 1 Online Subregion Updating

##  VI Hierarchical Planning Strategy

At this stage, it is imperative to consider organizing the access sequence of
these subregions to get an efficient global coverage path. Then, guided by the
global path, the visit order of local viewpoints is refined. As depicted in
Fig. 6, we employ a hierarchical planning paradigm to make a coarse-to-fine
decision on the next visit goal.

Figure 6: Illustration of hierarchical planning. The red arrows represent the
global coverage path $\mathcal{G}_{\text{global}}$ that sequentially enters
all subregions in active status, while the dashed lines indicate the optimized
local path $\mathcal{T}_{\text{local}}$ guided by the global path
$\mathcal{G}_{\text{global}}$. The local path visits all viewpoints within
$r_{\text{max}}$.

###  VI-A Global Coverage Path Planning

The objective of global coverage path planning is to find a global path
$\mathcal{G}_{\text{global}}$ that enters each subregion cell $S_{i}$ and
efficiently covers the entire unknown environment. Similar to [10], we also
formalize it as a variant of the TSP. Assuming there are $N_{\text{act}}$
subregions in active status, solving this TSP requires computing a
$(N_{\text{act}}+1)$-dimensional distance cost matrix
$\mathbf{M}_{\text{cost}}$, as the current robot pose $\mathbf{p}_{0}$ also
needs to be considered. Since the traveling cost
$\mathbf{L}_{\text{cost},i,j}$ of each subregion has already been computed in
Sec.V-C, $N_{\text{act}}\times N_{\text{act}}$ blocks can be directly filled
in:

$\mathbf{M}_{\text{cost}}(S_{i},S_{j})=\mathbf{L}_{\text{cost},i,j},i,j\in\left\\{1,2,\cdots,N_{\text{act}}\right\\}$
(4)

Note that $\mathbf{L}_{\text{cost},i,j}$ is computed based on the feasible
path between $\mathbf{C}_{\text{avg},i}$ of $S_{i}$ and
$\mathbf{C}_{\text{avg},j}$ of $S_{j}$, which is obtained in Sec.V-C by
executing the A* algorithm in the sparse graph $\mathcal{G}$. And only
$N_{\text{act}}$ additional calculations of the distance
$\mathbf{L}_{\text{cost},\mathbf{p}_{0},k}$ from the robot pose
$\mathbf{p}_{0}$ to each subregion cell are required:

$\mathbf{M}_{\text{cost}}(0,S_{k})=\mathbf{L}_{\text{cost},\mathbf{p}_{0},k},k\in\left\\{1,2,\cdots,N_{\text{act}}\right\\}$
(5)

Also adding a zero distance cost in the first column of
$\mathbf{M}_{\text{cost}}$ simplifies the TSP problem into an asymmetric TSP:

$\mathbf{M}_{\text{cost}}(S_{k},0)=0,k\in\left\\{1,2,\cdots,N_{\text{act}}\right\\}$
(6)

Thus far, the cost matrix $\mathbf{M}_{\text{cost}}$ is constructed, and the
TSP is then solved using a Lin-Kernighan-Helsgaun heuristic [25]. The output
of the solver is the sequence of accessing $N_{\text{act}}$ subregion cells
$\mathcal{G}^{*}_{\text{global}}=[S_{1},S_{2},\cdots,S_{N_{\text{act}}}]$,
where $\mathcal{G}^{*}_{\text{global}}$ denotes the optimal access sequence of
the subregion cells. This access sequence represents the shortest path that
can sequentially pass through all active subregions.

###  VI-B Local Path Optimization

The global coverage path identifies an efficient entry sequence for all active
subregion cells. However, there are still numerous viewpoints within each cell
that require visits to expand the frontiers. At this stage, our focus remains
on the cumulative traveling distance to visit these viewpoints, rather than
using a greedy strategy to evaluate their revenue.

To this end, we need to find the shortest visiting path
$\mathcal{T}^{*}_{\text{local}}$ that connects viewpoints while also
considering the endpoint constraints under the global coverage path
$\mathcal{G}^{*}_{\text{global}}$. Given the visiting order of subregion cells
$\mathcal{G}^{*}_{\text{global}}=[S_{1},S_{2},\cdots,S_{N_{\text{act}}}]$, the
local path is denote as:

$\mathcal{T}_{\text{local}}=[\mathbf{v}_{1}^{1},\mathbf{v}_{1}^{2},\cdots][\mathbf{v}_{2}^{1},\mathbf{v}_{2}^{2},\cdots]\ldots[\mathbf{v}^{1}_{N_{\text{act}}},\mathbf{v}^{2}_{N_{\text{act}}},\cdots]$
(7)

where $[\cdot]$ represents a local path segment in a subregion cell, and
$\mathbf{v}_{i}^{j}\in\mathbf{VP}_{i}$, which denotes the $j$-th viewpoint in
the $i$-th subregion. This formulation ensures that under the constraints of
$\mathcal{G}^{*}_{\text{\text{global}}}$, viewpoints in $S_{i}$ are visited
before those in $S_{i+1}$.

The above problem can be viewed as a TSP with endpoint constraints. Meanwhile,
guided by the global path, it is unnecessary to consider viewpoints that are
distant from the robot. Only viewpoints within a radius $r_{\text{max}}$ are
considered. Similar to the procedures outlined in Eq. 4 and Eq. 6, the A*
algorithm is used in $\mathcal{G}$ to search for feasible paths between
viewpoints to calculate distance cost $\mathbf{D}_{\text{cost}}$, and zero
costs are added to form the asymmetric TSP. The only difference is that in
addition to the distance cost $\mathbf{D}_{\text{cost},\mathbf{p}_{0}}$ from
the robot pose $\mathbf{p}_{0}$ to each viewpoint, a motion consistency cost
$c_{c}$ is also incorporated to discourage the back-and-forth behavior:

$\displaystyle\mathbf{M}_{\text{cost}}(0,\mathbf{v}_{l})=\mathbf{D}_{\text{cost},\mathbf{p}_{0},l}+\mu_{c}\cdot
c_{c}(\mathbf{p}_{0},\mathbf{v}_{l})$ , (8) $\displaystyle
l\in\left\\{1,2,\cdots,N_{\mathbf{v}}\right\\}$

$c_{c}(\mathbf{p}_{0},\mathbf{v}_{l})=\cos^{-1}\frac{(\mathbf{p}_{\mathbf{v}_{l}}-\mathbf{p}_{0})\cdot\mathbf{v}_{0}}{\|\mathbf{p}_{\mathbf{v}_{l}}-\mathbf{p}_{0}\|\|\mathbf{v}_{0}\|}$
(9)

where $\mathbf{v}_{l}$ represents the $l$-th viewpoint among a total of
$N_{\mathbf{v}}$ viewpoints in range $r_{\text{max}}$, $\mu_{c}$ is the tuning
factor, $\mathbf{v}_{0}$ is the current velocity of the robot, and
$\mathbf{p}_{\mathbf{v}_{l}}$ represents the position of the viewpoint
$\mathbf{v}_{l}$.

Note that the optimization problem mentioned above is not merely finding the
shortest path visiting each viewpoint. As stated in Eq. 7, it must satisfy the
endpoint constraint, meaning $\mathbf{v}_{i}^{j}$ must be visited before
$\mathbf{v}^{1}_{i+1}$. This constraint is included before solving with the
TSP solver.

After obtaining $\mathcal{T}^{*}_{\text{local}}$, the resulting discrete
waypoints consist of adjacent viewpoints. It is crucial to connect these
waypoints using the shortest passable path between viewpoints, which is then
transmitted to the motion planning module as the executable path to guide the
robot’s movement. The passable path between viewpoints has already been
determined during the calculation of $\mathbf{D}_{\text{cost}}$, ensuring that
this process incurs no additional computational overhead.

##  VII Experiments and Analysis

(a) Scene 1

(b) Scene 2

(c) Scene 3

Figure 7: The three simulation environments built in the Gazebo simulator.

###  VII-A Simulation Experiment

The simulation experiments are conducted on an Intel Core<EMAIL_ADDRESS>GHz
CPU with Ubuntu 20.04 LTS and 16 GB RAM. We evaluate the performance of our
method in simulation environments by comparing it with three state-of-the-art
methods, e.g., TARE[16], DSVP[5], and GBP [4]. The compared metrics include
exploration volume, traveling length, and algorithm runtime. Each method is
run 10 times in environments depicted in Fig. 7. All tests utilize the
exploration algorithm as the top-level decision module, while the middle
layers such as terrain analysis and obstacle avoidance employ an open-source
framework [26] from the Robotics Institute at Carnegie Mellon University. For
all approaches, the starting position is consistent, with a limited lidar
detection radius of 13 m and a maximum robot speed set to 2.0 m/s.

TABLE II: Simulation Results in Three Environments

Scene | Method |  Distance(m) |  Time Duration(s) |  Success Number
---|---|---|---|---
Avg | Std | Avg | Std
Scene 1 ($100\times 80$) | GBP |  $>$ 824.21 | - |  $>$ 1200 | - | 0
DSVP | 1440.85 | 67.18 | 866.27 | 30.27 | 10
TARE | 1348.60 | 52.11 | 822.78 | 33.41 | 10
Ours | 992.62 | 10.13 | 589.60 | 5.25 | 10
Scene 2 ($137\times 104$) | GBP |  $>$ 1235.55 | - |  $>$ 1600 | - | 0
DSVP | 1697.53 | 83.58 | 1057.24 | 61.34 | 10
TARE | 1568.47 | 60.38 | 885.91 | 29.93 | 10
Ours | 1439.82 | 29.39 | 819.56 | 13.22 | 10
Scene 3 ($80\times 80$) | GBP |  $>$ 606.27 | - |  $>$ 800 | - | 0
DSVP | 845.43 | 18.78 | 514.02 | 13.49 | 10
TARE | 778.55 | 46.53 | 492.99 | 28.46 | 10
Ours | 630.76 | 27.13 | 390.40 | 20.10 | 10

1) Exploration Performance: Table. II presents the specific experimental
results and statistics, while Fig. 8 visualizes the exploration progress
curves for each method. The exploration time limit for each scenario is set to
twice the maximum exploration time of our method. All methods except GBP cover
the entire exploration space. However, our proposed method consistently
outperforms the others in each environment. In the middle and late stages of
exploration, GBP struggles to relocate and switch exploration goals in open
areas, causing the robot to freeze in place and preventing it from fully
exploring the entire environment. DSVP employs a greedy exploration strategy,
selecting branches with the highest gain on the expanded tree. This often
leads to local optima and back-and-forth movements due to its focus on high-
gain regions while ignoring lower-gain areas.

Our approach, similar to TARE, addresses the TSP of arranging the visitation
order in the unknown region. However, we transform the coverage planning
problem into an open-loop asymmetric TSP, eliminating the need to return to
the starting point after each iteration. This approach avoids the necessity of
considering the return path from the last node to the start, which is more
typical in exploration tasks. In contrast, TARE addresses a closed-loop TSP,
where some regions are scheduled to be visited on the return trip. Since
exploration tasks generally do not require a return to the starting point, the
closed-loop TSP can result in an inefficient order of access, leading to
missed areas and subsequent backtracking. By employing an open-loop asymmetric
TSP, our strategy circumvents these issues, thereby enhancing exploration
efficiency. Moreover, in our method, neither the global coverage path nor the
local path accounts for information gain. We believe that introducing other
gains when constructing the cost matrix of the TSP may interfere with the
solving process and lead to suboptimal solutions. Since exploration is a long-
term task, minimizing exploration distance should be a primary goal.

This does not mean that information gain is trivial, but rather that once
sufficient information gain is achieved, its quantity can be disregarded. In
our method, subregion cells without viewpoints are excluded from global
coverage path planning because they lack sufficient information gain.
Similarly, all viewpoints are based on the collision-free sphere and can see
frontiers with ample information gain. From this perspective, our calculation
of information gain is implicitly included in the environment preprocessing
process.

(a) Scene 1

(b) Scene 2

(c) Scene 3

Figure 8: The exploration progress of four methods in each scene. The solid
line represents the average exploration progress, while the shading region in
color is formed by the upper-bound and lower-bound of the exploration
progress. The vertical dashed line represents the average exploration end
time.

(a) Scene 1

(b) Scene 2

(c) Scene 3

Figure 9: The algorithm runtime of three methods in each scene. The solid line
represents the average runtime. while the shading region in color is formed by
the upper-bound and lower-bound of the runtime.

2) Computational Efficiency: We also conduct a comprehensive analysis with
DSVP and TARE regarding algorithm runtime to highlight the advantages of our
proposed method, as shown in Fig. 9. Compared to TARE and DSVP, our method
exhibits the shortest runtime due to its rapid environmental preprocessing and
hierarchical planning strategy. DSVP suffers from long iteration times,
leading to delayed decision-making, particularly in large-scale environments.
The computational burden of DSVP primarily lies in evaluating the gain of
branches, which involves counting many unknown voxels and using ray tracing to
check sensor visibility. Consequently, more than 90% of the runtime is
consumed by this raycast-based gain computation. In contrast, our approach
samples viewpoints from nodes on the sparse graph using the collision-free
sphere, leveraging the number of frontier voxels observable from each
viewpoint as an information gain metric. This process requires minimal ray
tracing to assess visibility from the spherical region to nodes.

TARE addresses the access order of local planning horizons and global
subspaces. However, it generates a large number of viewpoints and subspaces
during exploration, especially in outdoor areas of Scene 1 and Scene 3. These
numerous viewpoints and subspaces need to be accessed sequentially by solving
the TSP, increasing the complexity of the solution. Our method continuously
decomposes the subregion to represent the unknown region with appropriate
resolution. All viewpoints are generated by collision-free spheres, minimizing
redundant viewpoints. Guided by the global-level coverage path, we solve the
visiting sequence for a small number of viewpoints, thereby reducing the
complexity of TSP solving.

###  VII-B Ablation Experiment

To study the specific impact of the greedy strategy on exploration efficiency,
we replace the hierarchical planning strategy with the mechanism of moving to
the nearest viewpoint as proposed in [6]. The results are shown in Table. III.
During exploration, the greedy strategy leads to frequent back-and-forth
movements between several viewpoints, reducing exploration efficiency.
Furthermore, this strategy can result in some areas being missed,
necessitating additional retracing to ensure complete coverage. Therefore, it
is crucial to assess the travel cost of visiting all areas to be explored.
This allows for arranging an efficient sequence for visiting different
regions, providing overall guidance for the entire exploration process.

We also investigate the impact of distance calculation methods on solving the
TSP. We use Euclidean distance to describe the distance between two areas
instead of the exact distance obtained by A* search. In constructing the TSP
distance matrix, we find that inaccurate distance measurements not only lead
to suboptimal TSP solutions but also cause frequent path switching, thereby
affecting overall optimization performance. Therefore, choosing an appropriate
distance measurement is crucial for effectively solving the TSP. Accurate
distance measurements can precisely reflect travel costs, ensuring a
reasonable visitation sequence.

TABLE III: Evaluation of Two Mechanisms on Exploration Performance

Scene | Method |  Distance(m) |  Time Duration(s)
---|---|---|---
Avg | Std | Avg | Std
Scene 1 | Greedy Strategy | 1167.10 | 15.12 | 763.30 | 32.20
Euclidean Distance | 1419.75 | 12.69 | 998.09 | 24.00
Scene 2 | Greedy Strategy | 1548.78 | 34.6 | 865.78 | 20.65
Euclidean Distance | 1606.65 | 27.5 | 949.68 | 24.98

###  VII-C Real-world Experiment

The real-world experiment is conducted in an underground garage with parked
vehicles, as shown in Fig. 1, with the size of the garage being 50 m $\times$
50 m. The environment includes moving vehicles and narrow passageways formed
between the parked vehicles. Fig. 1 also shows our real-world ground platform
with an AMD Ryzen 7<EMAIL_ADDRESS>GHz CPU and 32 GB RAM using Ubuntu 20.04 LTS.
The platform is topped with the RoboSense Hellos-16P LiDAR sensor. The Direct
LiDAR-Inertial Odometry [27] is used as the SLAM module. The max update range
radius of the UFOMap is set to 8 m, and the max speed of the robot is set to
1.0 m/s.

Fig. 1 shows the trajectory of the exploration result, and there are no
redundant paths during the whole exploration process. At the beginning of the
exploration, the environment is unknown to the robot. The robot explores
5673.25 m3 in 175.53 s and moves 157.53 m. Our method completely explores the
entire garage without leaving any areas. The experiment proves our method’s
efficiency and effectiveness.

##  VIII Conclusion

In this paper, we present an autonomous exploration system that is capable of
efficiently modeling unknown environments. The method introduces a rapid
environmental preprocess approach to provide basic information for subsequent
exploration planning, greatly reducing the computational load. We employ a
hierarchical planning strategy to arrange the entry sequence of each subregion
cell and determine the next viewpoint to be visited from coarse to fine,
striking a balance between exploration efficiency and computational cost. Our
method emphasizes the importance of global coverage path planning for long-
term exploration tasks rather than focusing solely on immediate gains,
especially in large-scale environments. This system is expected to be deployed
on wheeled-legged and aerial robots in the future, enabling autonomous
exploration in various complex terrains.

## References

  * [1] B. Yamauchi, “A frontier-based approach for autonomous exploration,” in _Proceedings 1997 IEEE International Symposium on Computational Intelligence in Robotics and Automation CIRA’97.’Towards New Computational Principles for Robotics and Automation’_. IEEE, 1997, pp. 146–151.
  * [2] A. Bircher, M. Kamel, K. Alexis, H. Oleynikova, and R. Siegwart, “Receding horizon” next-best-view” planner for 3d exploration,” in _2016 IEEE international conference on robotics and automation (ICRA)_. IEEE, 2016, pp. 1462–1468.
  * [3] T. Cieslewski, E. Kaufmann, and D. Scaramuzza, “Rapid exploration with multi-rotors: A frontier selection method for high speed flight,” in _2017 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_. IEEE, 2017, pp. 2135–2142.
  * [4] T. Dang, M. Tranzatto, S. Khattak, F. Mascarich, K. Alexis, and M. Hutter, “Graph-based subterranean exploration path planning using aerial and legged robots,” _Journal of Field Robotics_ , vol. 37, no. 8, pp. 1363–1388, 2020.
  * [5] H. Zhu, C. Cao, Y. Xia, S. Scherer, J. Zhang, and W. Wang, “Dsvp: Dual-stage viewpoint planner for rapid exploration by dynamic expansion,” in _2021 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_. IEEE, 2021, pp. 7623–7630.
  * [6] D. Duberg and P. Jensfelt, “Ufoexplorer: Fast and scalable sampling-based exploration with a graph-based planning structure,” _IEEE Robotics and Automation Letters_ , vol. 7, no. 2, pp. 2487–2494, 2022.
  * [7] Z. Meng, H. Qin, Z. Chen, X. Chen, H. Sun, F. Lin, and M. H. Ang, “A two-stage optimized next-view planning framework for 3-d unknown environment exploration, and structural reconstruction,” _IEEE Robotics and Automation Letters_ , vol. 2, no. 3, pp. 1680–1687, 2017.
  * [8] M. Kulich, J. Faigl, and L. Přeučil, “On distance utility in the exploration task,” in _2011 ieee international conference on robotics and automation_. IEEE, 2011, pp. 4455–4460.
  * [9] A. Dai, S. Papatheodorou, N. Funk, D. Tzoumanikas, and S. Leutenegger, “Fast frontier-based information-driven autonomous exploration with an mav,” in _2020 IEEE international conference on robotics and automation (ICRA)_. IEEE, 2020, pp. 9570–9576.
  * [10] B. Zhou, Y. Zhang, X. Chen, and S. Shen, “Fuel: Fast uav exploration using incremental frontier structure and hierarchical planning,” _IEEE Robotics and Automation Letters_ , vol. 6, no. 2, pp. 779–786, 2021.
  * [11] J. Huang, B. Zhou, Z. Fan, Y. Zhu, Y. Jie, L. Li, and H. Cheng, “Fael: fast autonomous exploration for large-scale environments with a mobile robot,” _IEEE Robotics and Automation Letters_ , vol. 8, no. 3, pp. 1667–1674, 2023.
  * [12] C. Connolly, “The determination of next best views,” in _Proceedings. 1985 IEEE international conference on robotics and automation_ , vol. 2. IEEE, 1985, pp. 432–435.
  * [13] Q. Bi, X. Zhang, J. Wen, Z. Pan, S. Zhang, R. Wang, and J. Yuan, “Cure: A hierarchical framework for multi-robot autonomous exploration inspired by centroids of unknown regions,” _IEEE Transactions on Automation Science and Engineering_ , 2023.
  * [14] B. Zhou, H. Xu, and S. Shen, “Racer: Rapid collaborative exploration with a decentralized multi-uav system,” _IEEE Transactions on Robotics_ , 2023.
  * [15] X. Zhao, C. Yu, E. Xu, and Y. Liu, “Tdle: 2-d lidar exploration with hierarchical planning using regional division,” in _2023 IEEE 19th International Conference on Automation Science and Engineering (CASE)_. IEEE, 2023, pp. 1–6.
  * [16] C. Cao, H. Zhu, H. Choset, and J. Zhang, “Tare: A hierarchical framework for efficiently exploring complex 3d environments.” in _Robotics: Science and Systems_ , vol. 5, 2021, p. 2.
  * [17] R. Reinhart, T. Dang, E. Hand, C. Papachristos, and K. Alexis, “Learning-based path planning for autonomous exploration of subterranean environments,” in _2020 IEEE International Conference on Robotics and Automation (ICRA)_. IEEE, 2020, pp. 1215–1221.
  * [18] F. Niroui, K. Zhang, Z. Kashino, and G. Nejat, “Deep reinforcement learning robot for search and rescue applications: Exploration in unknown cluttered environments,” _IEEE Robotics and Automation Letters_ , vol. 4, no. 2, pp. 610–617, 2019.
  * [19] Y. Cui, S. Ye, X. Xu, H. Sha, C. Wang, L. Lin, Z. Liu, R. Xiong, and Y. Wang, “Learning hierarchical graph-based policy for goal-reaching in unknown environments,” _IEEE Robotics and Automation Letters_ , 2024.
  * [20] D. Duberg and P. Jensfelt, “Ufomap: An efficient probabilistic 3d mapping framework that embraces the unknown,” _IEEE Robotics and Automation Letters_ , vol. 5, no. 4, pp. 6411–6418, 2020.
  * [21] H. Zhang, S. Wang, Y. Liu, P. Ji, R. Yu, and T. Chao, “Efp: Efficient frontier-based autonomous uav exploration strategy for unknown environments,” _IEEE Robotics and Automation Letters_ , 2024.
  * [22] N. Beckmann, H.-P. Kriegel, R. Schneider, and B. Seeger, “The r*-tree: An efficient and robust access method for points and rectangles,” in _Proceedings of the 1990 ACM SIGMOD international conference on Management of data_ , 1990, pp. 322–331.
  * [23] Y. Ren, F. Zhu, W. Liu, Z. Wang, Y. Lin, F. Gao, and F. Zhang, “Bubble planner: Planning high-speed smooth quadrotor trajectories using receding corridors,” in _2022 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_. IEEE, 2022, pp. 6332–6339.
  * [24] B. Tang, Y. Ren, F. Zhu, R. He, S. Liang, F. Kong, and F. Zhang, “Bubble explorer: Fast uav exploration in large-scale and cluttered 3d-environments using occlusion-free spheres,” in _2023 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_. IEEE, 2023, pp. 1118–1125.
  * [25] K. Helsgaun, “An effective implementation of the lin–kernighan traveling salesman heuristic,” _European journal of operational research_ , vol. 126, no. 1, pp. 106–130, 2000.
  * [26] C. Cao, H. Zhu, F. Yang, Y. Xia, H. Choset, J. Oh, and J. Zhang, “Autonomous exploration development environment and the planning algorithms,” in _2022 International Conference on Robotics and Automation (ICRA)_. IEEE, 2022, pp. 8921–8928.
  * [27] K. Chen, R. Nemiroff, and B. T. Lopez, “Direct lidar-inertial odometry: Lightweight lio with continuous-time motion correction,” _2023 IEEE International Conference on Robotics and Automation (ICRA)_ , pp. 3983–3989, 2023.


# Gamma-ray observation towards the young massive star cluster NGC 6618 in the
M17 region

Bing Liu,1,2,3 Rui-zhi Yang,1,2,3 Zhiwei Chen4
1Department of Astronomy, School of Physical Sciences, University of Science
and Technology of China, Hefei, Anhui 230026, China
2CAS Key Labrotory for Research in Galaxies and Cosmology, University of
Science and Technology of China, Hefei, Anhui 230026, China
3School of Astronomy and Space Science, University of Science and Technology
of China, Hefei, Anhui 230026, China
4Purple Mountain Observatory, Chinese Academy of Sciences, Yuanhua Road 10,
210023 Nanjing, China E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

Young massive clusters have been established as a new population of gamma-ray
sources and potential cosmic ray (CR) accelerators. In this paper, we report
the detection of gamma-ray emissions near the young star cluster NGC 6618,
which is one of the youngest star clusters in our Galaxy. The detected gamma-
ray emissions can be divided into two components. One component is point-like
and reveals harder spectrum, while the other is extended and with softer
spectrum. Such spectral features are significantly different from other young
massive clusters and may be due to the propagation effects of CRs accelerated
in NGC 6618.

###### keywords:

cosmic rays – gamma-rays: ISM – open clusters and associations: individual:
M17 (NGC 6618)

††pubyear: 2021††pagerange: Gamma-ray observation towards the young massive
star cluster NGC 6618 in the M17 region–Gamma-ray observation towards the
young massive star cluster NGC 6618 in the M17 region

##  1 Introduction

Young massive star clusters (YMCs) are established as a new population of
gamma-ray sources in recent years and are believed to be able to account for
the acceleration of a significant part of Galactic cosmic rays (CRs)
(Aharonian et al., 2019). Several such systems are detected in GeV or TeV
gamma-ray band, with significant spatial extensions up to more than 50 pc and
a universal hard gamma-ray spectrum which can be described by a powerlaw
function with an index of $2.2-2.4$. More interestingly, the derived CR
spatial distribution reveals a $\frac{1}{r}$ profile, which is consistent with
a continuous injection of CRs. Considering the size of the source, the CRs are
believed to be injected in the time scale of $\ga 10^{5}$ years, which is
longer than the acceleration phase of supernova remnants (SNRs). Moreover, the
massive star winds, which can be as high as $\ga 3000~{}\rm km\,s^{-1}$ in
O-type stars and Wolf-Rayet stars, are powerful enough to accelerate CRs to
high energies nearly in the whole lifetime of the massive stars (Cesarsky &
Montmerle, 1983). In this regard, the CRs in YMCs are more likely accelerated
by "live" stars rather than the "dead" stars. In an observational point of
view, several YMCs have been detected in gamma-rays, such as Cygnus
Cocoon(Ackermann et al., 2011), NGC 3603 (Yang & Aharonian, 2017; Saha et al.,
2020), Westerlund 1 (Abramowski et al., 2012), Westerlund 2 (Yang et al.,
2018), RSGC 1(Katsuta et al., 2017; Sun et al., 2020a), W40 (Sun et al.,
2020b) and W43 (Yang & Wang, 2020).

The YMC NGC 6618, with an age less than 3 Myr (Jiang et al., 2002), is one of
the youngest massive star clusters in our Galaxy, and contains hundreds of
stars earlier than B9 including dozens of O stars (Chini et al., 1980;
Hoffmeister et al., 2008). The OB stars in the cluster, especially the most
massive ones such as the O4 binary system (CEN 1a + CEN 1b; Hoffmeister et
al., 2008), provide the ionizing power for the H ii region, M17 nebula (also
known as Omega Nebula or W38), which is located in the Sagittarius spiral arm
(Elmegreen et al., 1979; Reid et al., 2019) at a distance $\approx$ 2.0 kpc
(Xu et al., 2011; Chibueze et al., 2016). Moreover, the cluster is also
thought to be the source of energetic feedback to the dense molecular clouds
(MCs) surrounding M17 nebula (Nguyen-Luong et al., 2020).

The MC associated with M17 was first detected in CO molecular line emission at
velocity around $20\,\mathrm{km\,s^{-1}}$, and the line emission extends along
the Galactic plane for at least 85 pc and has a total mass in the order of
$10^{6}\,{\rm M}_{\sun}$ (Elmegreen & Lada, 1976), comparable to the typical
definition of the giant MC (GMC). The M17 GMC is usually split into three
regions according to the relative location with respect to NGC 6618. The cloud
in the north, i.e., M17 North cloud, has the smallest mass and size (Wilson et
al., 2003). The adjacent cloud in the southwest, i.e., M17 SW cloud, is the
densest region and highly clumpy (Stutzki & Guesten, 1990; Nguyen-Luong et
al., 2020), and contributes most of the forming stars of low- to high-mass
(Povich et al., 2009; Chen et al., 2012; Chen et al., 2013, 2015, 2021;
Ramírez-Tannus et al., 2017; Lim et al., 2020). The cloud extending to the
very southwest, i.e., M17 SWex, appears as an infrared dark cloud and hosts
sites of low- to intermediate-mass star formation (Povich & Whitney, 2010;
Povich et al., 2016; Yin et al., 2022). By comparing the dense gas properties
of M17 SW and M17 SWex, Nguyen-Luong et al. (2020) found that the clumps in
M17 SW are denser, more compact, and more gravitationally bound than those in
M17 SWex.

If the massive stars in NGC 6618 can accelerate protons to very high energies,
then these protons could illuminate the clouds around them via proton-proton
inelastic collisions. Recently released Fermi-LAT 10-year Source Catalog
(4FGL-DR2, Ballet et al., 2020) revealed a gamma-ray point source that located
in the direction of M17 region (named as 4FGL J1820.4-1609c), which is
possibly associated with NGC 6618. In such a complex region as described
above, the origin of these gamma-rays is still an intriguing puzzle and
requires a careful and comprehensive investigation. In the following section,
we present the details of the analysis of Fermi-LAT data and the results. Then
in Sect. 3, we investigate the gas distributions around M17. Next in Sect.4,
we combine the multiwavelength observations with gamma-ray results to explore
the radiation mechanism of the gamma-rays. Finally, we summarise the main
conclusions and implications of our study in Sect.5.

##  2 Fermi-LAT data analysis and results

To study the gamma-ray emission in M17 region, we collected about 12.5 years
(from 2008-08-04 15:43:36 (UTC) to 2021-04-16 23:02:55 (UTC)) of Fermi-LAT
Pass 8 data, and used Fermitools from Conda
distribution111https://github.com/fermi-lat/Fermitools-conda/ and applied the
latest version of the instrument response functions P8R3_SOURCE_V3. Given the
energy-dependent point-spread function (PSF) of the LAT data, we used
different datasets to optimize the spatial and spectral analyses. The region
of interest (ROI), centered at M17 (R.A.$=275.195^{\circ}$,
Dec.$=-16.172^{\circ}$, J2000) is adjusted for each dataset correspondingly.
The source models, generated by
make4FGLxml.py222https://fermi.gsfc.nasa.gov/ssc/data/analysis/user/make4FGLxml.py,
consist of the sources in 4FGL-DR2 within the ROI enlarged by $10^{\circ}$,
the Galactic diffuse background emission (_gll_iem_v07.fits_), and isotropic
emission (_iso_P8R3_SOURCE_V3_v1.txt_).

###  2.1 Spatial analysis and results

In order to study possible energy dependent morphology in M17 region, we
conduct separate spatial analyses for gamma-rays within different energy
ranges, including 0.5$-$5 GeV, 5$-$500 GeV, and 0.5$-$500 GeV, hereinafter
referred as LE, HE and ALL datasets, respectively. Specific parameters of the
cuts and details of each dataset are provided in Table.1.

Following procedures were applied to the LE and HE datasets respectively, in
which the spectral type of the tested sources are all set to be simple
PowerLaw (PL). At first, we used gtlike to optimize the source models. Since
4FGL J1820.4-1609c is very likely associated with NGC 6618 and cannot be
treated as background. Thus we removed 4FGL J1820.4-1609c from the optimized
background source models, and generated the residual test statistic (TS) maps.
As shown in Fig.1, the residual gamma-rays of LE dataset are very diffusive
and shifted to the west of M17, meanwhile, the gamma-rays of HE dataset are
rather point-like and in the direction of NGC 6618. Therefore, to find out
whether the residual emission is extended or not, we added a point-like source
or an extended source (using RadialGaussian model) to the modified source
model, from which 4FGL J1820.4-1609c are deleted. The added point or the
center of the gaussian disk is set at the peak position of the residual TS map
for each dataset, which is presented in Table 2. The radius ($\sigma_{\rm
disk}$) of the disk varies from $0.1^{\circ}$ to $0.5^{\circ}$ with a step of
$0.05^{\circ}$. Then we fitted above modified models to calculate the
corresponding $\rm TS_{\rm ext}$, defined as $-2\log({\mathcal{L}}_{\rm
p}/{\mathcal{L}}_{\rm ext})$. The $L_{\rm p}$ is the maximum likelihood of the
model with a point source added at the peak position, and the
$\mathcal{L}_{\rm ext}$ is the acquired maximum likelihood of the tested
gaussian disk models. For LE dataset, the maximum $\rm TS_{\rm ext}$ (290.6)
is obtained when $\sigma_{\rm disk}$=$0.35^{\circ}$, meaning the significance
of the extension is about 17.0. However, for HE dataset, gaussian disk models
show no improvement compared to the point source assumption. The specific
results for the spatial analyses of LE and HE datasets are listed in Table 2.

Finally, we analysed ALL dataset in order to test whether one source with the
spatial model obtained from LE or HE dataset (1-source hypothesis) or two
sources with the spatial models acquired from both LE and HE datasets
(2-source hypothesis) could represents the overall data better. By comparing
the maximum log(likelihood) of these models, we obtained the significance of
2-source hypothesis versus 1-source hypotheses or the non-modified 4FGL-DR2
model is $\sim 9\sigma$.

###  2.2 Spectral analysis and results

From above analysis, we noticed that the morphology of gamma-ray emission
around M17 region is energy-dependent and these gamma-rays are better
represented as two separate components. In the following spectral analysis
applied to ALL dataset, we replace 4FGL J1820.4-1609c with two sources: one
source is an extended source with $\sigma_{\rm disk}=0.35^{\circ}$ lying to
the west of NGC 6618 (hereafter referred as Src A), and the other is a point-
like source that is coincided with NGC 6618 (hereafter referred as Src B).

To find out the spectral shape of Src A and Src B, we performed likelihood-
ratio test for spectral models including LogParabola (LogP), PLSuperExpCutoff2
(PLEC), BrokenPowerLaw (BPL), in which the PL model is the null hypothesis.
The formulae and free parameters of these spectral model are presented in
Table.3. The significance of the tested model $\sigma_{\rm model}$ is defined
as $\sqrt{-2\log({\mathcal{L}}_{\rm PL}/{\mathcal{L}}_{\rm model})}$. Firstly,
we changed the spectral type of Src A from PL to LogP, PLEC and BKL
respectively, keeping the spectral model of Src B to be PL. Next, we set the
spectral model of Src A to be PL, and change the spectral type of Src B from
PL to LogP, PLEC, and BPL, respectively. As shown in Table.4, those more
complicated spectral models do not improve the overall fitting results
($\sigma_{\rm model}<3$), and the simple PL model is capable to represent
their spectral shape. The best-fit results of the PL model for Src A and Src B
are listed in Table.5.

The spectral energy distributions (SEDs) of Src A and Src B as shown in Fig.2
were extracted from the maximum likelihood analysis of source class events in
nine logarithmically spaced energy bins within 0.5–500 GeV. Expect the energy
range of the selected events, other parameters applied for the data
preparation are the same as the ALL datasets. During the fitting process, the
free parameters only include the normalization parameters of the sources with
the significance $\geq 5\sigma$ within $10^{\circ}$ from ROI center as well as
the Galactic and isotropic diffuse background components, while all the other
parameters are fixed to their best-fit values from above analysis of ALL
dataset in which the spectral shape of Src A and Src B are assumed to be PL.
In addition to the uncertainties caused by the statistics and the effective
area (green error bars in Fig.2), we also estimated the uncertainties caused
by the imperfection model of the Galactic diffuse background (blue error bars
in Fig.2) by artificially varying its normalization by $\pm 6\%$ from the
best-fit value of each energy bin and recording flux deviation of the source
due to above changes as the systematic error, following the method from Abdo
et al. (2009).We note that in lower energy range (0.5–1 GeV, first bin of the
SEDs), the fluxes of both Src A and Src B drop dramatically when the Galactic
diffuse background was artificially enhanced.

Table 1: Description of the datasets for spatial analysis

Name | ALL | LE | HE
---|---|---|---
Emin (GeV) | 0.5 | 0.5 | 5
Emax(GeV) | 500 | 5 | 500
ROI | $20^{\circ}\times 20^{\circ}$ | $20^{\circ}\times 20^{\circ}$ | $10^{\circ}\times 10^{\circ}$
zmax | $90^{\circ}$ | $90^{\circ}$ | $105^{\circ}$
evclass | 128 | 128 | 128
evtype | 16/32 | 32 | 16
Nebins | 30 | 10 | 16
DATA_QUAL | $>0$ | $>0$ | $>0$
LAT_CONFIG | 1 | 1 | 1
freea/allb |  $10^{\circ}$/ $25^{\circ}$ |  $10^{\circ}$/ $25^{\circ}$ |  $5^{\circ}$/ $15^{\circ}$

a The radius in degrees from the center of ROI within which the spectral
parameters are free to change for sources with significance $>5$.
b The radius in degrees from the center of ROI within which sources are
included in the XML file.

Table 2: Spatial analysis results of the LE and HE datasets

Name | LE | HE
---|---|---
Position (Ra, Dec) | (274.65∘ ,-16.30∘) | (275.120∘ ,-16.20∘)
$\mathcal{L}_{\rm 0}$ | 645931.3 | -36584.1
$\mathcal{L}_{\rm p}$ | 646128.2 | -36579.1
$\mathcal{L}_{\rm ext}$ | 646273.5 | -36578.8
$\sigma_{\rm disk}$ (∘ ) | $0.35\pm 0.05$ | $0.10\pm 0.05$
${\rm TS}_{\rm ext}$ | 290.6 | 0.3

The position refers to the peak location of the residual TS map for each
dataset. $\sigma_{\rm disk}$ is the radius of gaussian disk of the maximum
-log(likelihood) among the extension test models in which $\sigma_{\rm disk}$
varies from $0.1^{\circ}$ to $0.5^{\circ}$ with a step of $0.05^{\circ}$ and
the center of the disk located at the corresponding peak position. The
$\mathcal{L}_{\rm 0}$ is the likelihood of the non-modified 4FGL-DR2 model.
The $\mathcal{L}_{\rm p}$ is the likelihood of the model with 4FGL
J1820.4-1609c removed and a point source added at the peak position. The
$\mathcal{L}_{\rm ext}$ is the -log(likelihood) of best-fit gaussian disk
model. ${\rm TS}_{\rm ext}$ =$-2\log({\mathcal{L}}_{\rm p}/{\mathcal{L}}_{\rm
ext})$.

Figure 1: Smoothed residual TS maps centered at M17 for LE dataset (left) and HE dataset (right), respectively. The green crosses show the position of nearby sources listed in 4FGL-DR2, and the magenta diamond represents 4FGL J1820.4-1609c, which is removed from the fitted background model. The white dashed circle and crosses illustrate the nearby supernova remnants and pulsars, respectively. The black dashed circle illustrate the best-fit gaussian disk model (Src A) from fitting LE dataset, and the black cross represents the position of the point source (Src B) for fitting HE dataset. The maps are overlaid with contours generated from 13CO ($J$=1-0) integrated intensity map in the velocity range 10-34 $\mathrm{km\,s^{-1}}$ at ten linear scale levels between 15 and 150 $\mathrm{K\,km\,s^{-1}}$. See Sect. 2.1 for details.  Table 3: Formulae for gamma-ray spectral analysis Name | Formula | Free parameters
---|---|---
PL |  $\frac{\mathrm{d}N}{\mathrm{d}E}$ = $N_{0}{\left(\frac{E}{E_{0}}\right)}^{-\Gamma}$ |  $N_{0}$, $\Gamma$
LogP |  $\frac{\mathrm{d}N}{\mathrm{d}E}$ = $N_{0}\left(\frac{E}{E_{\mathrm{b}}}\right)^{-\Gamma-\beta\log(\left(\frac{E}{E_{\mathrm{b}}}\right)}$ |  $N_{0}$, $\Gamma$, $\beta$
PLEC |  $\frac{\mathrm{d}N}{\mathrm{d}E}$ = $N_{0}\left(\frac{E}{E_{0}}\right)^{-\Gamma_{1}}\exp\left({-\left({\frac{E}{E_{\mathrm{cut}}}}\right)^{\rm b}}\right)$ |  $N_{0}$, $\Gamma$, $E_{\mathrm{cut}}$, b
BPL |  $\frac{\mathrm{d}N}{\mathrm{d}E}$ = $\begin{cases}N_{0}{\left(\frac{E}{E_{\mathrm{b}}}\right)}^{-\Gamma_{1}}&\mbox{: }E<E_{\mathrm{b}}\\\ N_{0}\left(\frac{E}{E_{\mathrm{b}}}\right)^{-\Gamma_{2}}&\mbox{: }E>E_{\mathrm{b}}\end{cases}$ |  $N_{0}$, $\Gamma_{1}$, $\Gamma_{2}$, $E_{\rm b}$
Table 4:  The significance of the tested spectral models ($\sigma_{\rm
model}$) for Src A and Src B

Spectral Model | LogP | PLEC | BPL
---|---|---|---
Src Aa | 2.3 | 2.3 | 2.2
Src Bb | 2.5 | 2.3 | 2.3

a The spectral type of Src B is PL
b The spectral type of Src A is PL

Table 5: The best-fit results of the PL model for Src A and Src B

| $\Gamma$ | Flux | Luminosity a | significance b
---|---|---|---|---
|  | ($\,\rm{ph}\,{\rm cm}^{-2}\,{\rm s}^{-1}$) | ($\rm{erg}\,{\rm s}^{-1}$) | $\sigma$
Src A | $2.73\pm 0.08$ | $1.58\times 10^{-8}$ | $1.14\times 10^{33}$ | 15
Src B | $2.41\pm 0.06$ | $5.67\times 10^{-9}$ | $5.67\times 10^{32}$ | 12

a The distance of Src A and Src B is assumed to be 2 kpc.
b The significance of each source $\simeq\sqrt{\rm TS}$.

Figure 2:  SEDs of Src A (left) and Src B (right) respectively. For each data
point, the green error bar indicates the uncertainty caused by the statistics
and the effective area, and the blue error bar shows the uncertainty caused by
the imperfection model of the Galactic diffuse background. The best-fit PP
models (black lines) and IC models (red lines) for Src A and Src B are
represented. See Sect. 2.2 and Sect. 4 for details.

##  3 Gas distributions near M17

In this work, we analysed the large scale 12CO/13CO ($J$=1-0) molecular line
data toward the M17 GMC, as a part of the Milky Way Imaging Scroll Painting
(MWISP) project. The details of MWISP project are referred to Su et al.
(2019). The 12CO/13CO ($J$=1-0) molecular data used in this work span from the
galactic longitude $13\degr$ to $16\degr$ and the galactic latitude
$-1\aas@@fstack{\circ}5$ to $0\aas@@fstack{\circ}5$. After fitting the
baseline and calibrating the main beam efficiency, the reduced 3D data cubes
with a grid spacing of $30\arcsec$ have a typical root mean square noise level
of $\sim 0.5\,\mathrm{K}$ for 12CO ($J$=1-0) transition line and $\sim
0.3\,\mathrm{K}$ for 13CO ($J$=1-0) transition line at a channel width of
$0.16\,\mathrm{km\,s^{-1}}$. The spatial resolution of the CO data is $\sim
50\arcsec$.

The large-scale 12CO and 13CO gas distributions over the area from the
galactic longitude $13\aas@@fstack{\circ}67$ to $15\aas@@fstack{\circ}39$ and
the galactic latitude $-0\aas@@fstack{\circ}87$ to $-0\aas@@fstack{\circ}27$
were studied earlier in Nguyen-Luong et al. (2020); however, the spatial
extent of Src A is not fully covered by this research. The bulk of the M17 GMC
emission is seen in the velocity range $10-30\,\mathrm{km\,s^{-1}}$ from the
Sagittarius arm, while the components at $\sim 37-38\,\mathrm{km\,s^{-1}}$ and
$\sim 57\,\mathrm{km\,s^{-1}}$ are from the spiral arms at larger distances
(Zhan et al., 2016; Nguyen-Luong et al., 2020). Therefore, we only considered
the emission around $10-30\,\mathrm{km\,s^{-1}}$ for the MWISP CO data.

Figure 3: The integrated intensity map of 13CO ($J$=1-0) line emission in the
velocity range 25-33 $\mathrm{km\,s^{-1}}$ (left) and 17-25
$\mathrm{km\,s^{-1}}$ (right) respectively. The white circle represents the
extended gamma-ray component (Src A), and the white cross shows the location
of point-like component (Src B). See Sect. 3 for details.

As shown in Fig.3, Src B is much closer to the M17 SW cloud (hereafter
referred to as cloud B) than to the M17 north cloud. Cloud B shows a peak
intensity at the velocity $V_{\mathrm{LSR}}\approx 20\,\mathrm{km\,s^{-1}}$,
and a velocity range from 10 to $28\,\mathrm{km\,s^{-1}}$ (see Fig.4). At
velocity greater than $28\,\mathrm{km\,s^{-1}}$, the CO ($J$=1-0) line
emission shows an extended distribution at a level of $T_{\mathrm{MB}}\sim
2-3\,\mathrm{K}$. . The molecular gas within the region of Src A (hereafter
referred to as cloud A, the white circle in Fig.3) shows velocity in the range
from 10 to $34\,\mathrm{km\,s^{-1}}$, much broader than the velocity range of
cloud B (see Fig.4). For cloud A, the peak intensity of the 13CO ($J$=1-0)
line emission is at $V_{\mathrm{LSR}}\approx 27\,\mathrm{km\,s^{-1}}$.
However, both cloud B and M17 SWex cloud show peak intensity at
$V_{\mathrm{LSR}}\approx 20\,\mathrm{km\,s^{-1}}$. The large difference in
$V_{\mathrm{LSR}}$, together with the far distance from Src A to cloud B and
M17 SWex cloud, indicate that Src A is not associated with the molecular gas
at the velocity $V_{\mathrm{LSR}}\approx 20\,\mathrm{km\,s^{-1}}$. Thus, we
suggest that Src A is associated with cloud A with a peak intensity at
$V_{\mathrm{LSR}}\approx 27\,\mathrm{km\,s^{-1}}$.

Figure 4: Mean CO spectra for molecular gases (clouds A and B) that are
possibly associated with Src A (top) and Src B (bottom).

The bulk masses of clouds A and B depend strongly on the adopted velocity
range of the 12CO line emission. For cloud B that possibly associated with Src
B, a velocity range $13-25\,\mathrm{km\,s^{-1}}$ is assumed for the 12CO line
emission. A wider range in velocity would definitely include more molecular
gas. The molecular gas mass estimated from the velocity range
$13-25\,\mathrm{km\,s^{-1}}$ represents the lower limit of the mass of cloud
B. We see from Figure 4 that the 12CO line spectrum averaged over Src A shows
line width of $18\,\mathrm{km\,s^{-1}}$ and equal contributions from clouds A
and B. Since it is impossible for us to confirm an exact range in velocity for
cloud A that is likely related to Src A, here we adopt a velocity range
$25-34\,\mathrm{km\,s^{-1}}$ for cloud A. This velocity range would miss a
fraction of molecular gas with velocity exceeding this range, but preclude
mostly the molecular gas at velocity about $40\,\mathrm{km\,s^{-1}}$ from the
Scutum arm. Therefore, the integrated intensity of 12CO line emission ($W_{\rm
CO}$) over the ranges $25-34\,\mathrm{km\,s^{-1}}$ and
$13-25\,\mathrm{km\,s^{-1}}$ are obtained for clouds A and B, respectively.

Next, we used the nominal $X$-factor, $2\times
10^{20}\,\mathrm{cm^{-2}}$/($\mathrm{K\,km\,s^{-1}}$), for the conversion from
$W_{\mathrm{CO}}$ to the molecular hydrogen column density
$N(\mathrm{H_{2}})$. For cloud A, 13CO line emission shows a prominent cavity
of little emission toward the center of cloud A, while 12CO line emission
toward the same cavity is still significant. Thus, for cloud A, the average
$N(\mathrm{H_{2}})$ estimated from the 13CO line emission (hereafter referred
as $N(\mathrm{H_{2}})_{\rm{}^{13}CO}$) is $8.1\times
10^{21}\,\mathrm{cm^{-2}}$, slightly lower than the average
$N(\mathrm{H_{2}})$ estimated from 12CO line emission
($N(\mathrm{H_{2}})_{\rm{}^{12}CO}$), which is $1.0\times
10^{22}\,\mathrm{cm^{-2}}$. As for cloud B, the
$N(\mathrm{H_{2}})_{\rm{}^{12}CO}$ is $4.9\times 10^{22}\,\mathrm{cm^{-2}}$,
which should be severely underestimated because the 12CO line emission is
optically thick. In comparison, the estimated
$N(\mathrm{H_{2}})_{\rm{}^{13}CO}$ of cloud B can reach to $3.4\times
10^{23}\,\mathrm{cm^{-2}}$ 333This method makes use of both 12CO and 13CO line
emissions and assumes the local thermal equilibrium state for the two
molecules. Detailed description of this method is referred to Chen et al.
(2022) and references therein..

Assuming a distance of 2.0 kpc, the radius of cloud B is $\sim 2.5$ pc. Thus,
the total gas mass of cloud B is $2.2\times 10^{4}\,{\rm M}_{\sun}$ according
to its $N(\mathrm{H_{2}})_{\rm{}^{12}CO}$ or $1.5\times 10^{5}\,{\rm
M}_{\sun}$ according to its $N(\mathrm{H_{2}})_{\rm{}^{13}CO}$. Nguyen-Luong
et al. (2020) estimated a total gas mass of $1.43\times 10^{5}\,{\rm
M}_{\sun}$ for the molecular gas near to the YMC NGC 6618, including cloud B
and the less massive M17 North cloud. Therefore, we adopt a total gas mass of
$1.5\times 10^{5}\,{\rm M}_{\sun}$. Then the corresponding hydrogen number
density $n(\mathrm{H})$ of cloud B is $4.4\times 10^{4}\,\mathrm{cm^{-3}}$,
assuming a thickness of 5 pc. Similarly, the total gas mass of cloud A is
about $1.3\times 10^{5}\,{\rm M}_{\sun}$, comparable to that of cloud B.
Assuming a thickness of 5 pc for cloud A, a value as the same as cloud B, the
hydrogen number density $n(\mathrm{H})$ of cloud A is estimated to be
$1.3\times 10^{3}\,\mathrm{cm^{-3}}$. However, if the thickness of cloud A is
comparable to the projected size on the plane of sky, i.e., be as large as 25
pc, it will yield a lower density, $\sim 2.6\times 10^{2}\,\mathrm{cm^{-3}}$.
Thus, we suggest that $n(\mathrm{H})$ of cloud A is likely between $\sim
2.6\times 10^{2}\,\mathrm{cm^{-3}}$ and $\sim 1.3\times
10^{3}\,\mathrm{cm^{-3}}$.

In short, the molecular gas components possibly associated with Src A and Src
B, clouds A and B, are very different in density, size, and kinematics. At a
$V_{\mathrm{LSR}}\sim 27\,\mathrm{km\,s^{-1}}$, cloud A is a very extended
molecular cloud of moderate density. Cloud B is a compact dense molecular
cloud with a $V_{\mathrm{LSR}}\sim 20\,\mathrm{km\,s^{-1}}$.

##  4 Discussion

To investigate the possible radiation mechanisms of the gamma-rays in M17
region, we fit the SEDs of Src A and Src B with both leptonic scenario, i.e.,
the inverse Compton scattering (hereafter referred as IC model) and hadronic
scenario, i.e., proton-proton inelastic collision (hereafter referred as PP
model). The fitting was performed using Naima
package444http://naima.readthedocs.org/en/latest/index.html (Zabalza, 2015),
which includes tools to perform Markov Chain Monte Carlo (MCMC) fitting of
nonthermal radiative processes to the data and allows us to implement
different functions. Here the distribution function of the parent particles
was assumed to be simple power law. The best-fit results are illustrated in
Fig.2, in which the black lines are the results of PP model fitting and the
red lines are results of IC model fitting. The corresponding parameters are
listed in Table.6.

Around M17 region, there are several known pulsars and SNRs (see Fig.1), but
none of them overlaps with Src A or Src B. Although we cannot formally rule
out the possibility that the emissions are related with some unknown pulsars
or SNRs, the extension of Src A and the hard spectrum of Src B can exclude the
possibilities that those sources are pulsars. However, the pulsar wind nebula
related with the unknown pulsars can be a possible explanation. In addition,
considering the spatial correlation of both Src A and Src B with dense gas, it
is possible that the CRs escaped from SNRs in the vicinity interacting with
dense gas produced the detected gamma-rays.

Another scenario is associating Src A and Src B with the YMC NGC 6618. NGC
6618 contains more than 40 OB stars (Hoffmeister et al., 2008), the total wind
power of these stars can be estimated using the estimations in Domingo-
Santamaría & Torres (2006), which amounts to more than $5\times 10^{37}{\rm
erg\,s^{-1}}$. Taking into account the age of 500 kyr, the total energy
injected by NGC 6618 is about $10^{51}$ erg, which is comparable to a typical
supernova explosion. Thus, this system is powerful enough to accelerate enough
CRs to account for the detected gamma-ray emissions. Furthermore, the hard
spectrum of Src B is also similar to other YMC systems. The soft spectrum of
Src A is unique compared with other YMCs. Generally speaking, if we assume the
CRs are injected continuously by the YMC, the CR energy spectrum in the
vicinity can be expressed as
$F\left(E\right)\sim\frac{Q\left(E\right)}{D\left(E\right)}$, where
$Q\left(E\right)$ and $D\left(E\right)$ is the injection spectrum and
diffusion coefficient, respectively (Aharonian & Atoyan, 1996). Therefore, a
natural explanation for the different spectral indices of Src A and Src B is a
different $D\left(E\right)$ in the corresponding regions. Indeed, due to the
higher density and stronger magnetic field inside the YMCs, the environment
can be significantly different from the interstellar medium (ISM) in the
Galaxy. The MHD turbulent cascade in the ISM can be damped out effectively and
CRs will stream along field lines and transport via a process of field line
random walk. Such process have been studied in detail in starburst galaxies in
Krumholz et al. (2020), and the environments in YMCs are comparable to these
starburst galaxies. In this case, the effective diffusion coefficient is
energy independent at lower energy (below some critical energy $E_{\rm br}$),
thus the propagated CR spectra are the same as the injected spectra. This
provide an natural explanation for hard spectrum in Src B and in other YMC
systems. If we assume a continuous injection, the CR energy density scales as
$\frac{1}{r}$, where $r$ is the distance to the CR source. Thus, the CR
density can be much smaller in Src A than in Src B, if we assume the CRs are
injected from NGC 6618. For Src A, the CR transport may be still dominated by
the CR scattering with the MHD turbulence in the ISM. In this case,
$D\left(E\right)$ scales as $E^{0.33...0.5}$, which predicts a softer
spectrum. The difference of indices of Src A and Src B is about 0.3, which is
also consistent with the index of $D\left(E\right)$ in the ISM.

In such a scenario all the gamma-ray emissions surrounding YMCs should contain
two regions, in one of which the CR transport is dominated by streaming along
field line and reveal a harder spectrum (inner region); in another the CR
transport is dominated by scattering with MHD turbulence in the ISM (outer
region). For other YMC systems, the "outer region" has not been detected yet,
one reason may be that for those more powerful systems such as Cygnus cocoon,
Westerlund 2, and NGC 3603, due to the larger wind power and longer age, the
CR density therein is significantly higher, thus the "inner region" in which
CRs transport by streaming along field lines is much larger than that in NGC
6618, and occupy the whole dense regions near the YMCs, thus dominate the
produced gamma-ray emissions. In weaker systems such as W40 and W43, the
"outer region" in which CR transport dominated by diffusion in MHD turbulence
in ISM are too weak to be distinguished from the diffuse gamma-ray background.

In the case of "inner region", the CR streaming velocity is dominated by the
balance between the damping rate and the CR streaming instability, which
results a energy dependent diffusion at higher energy and predict a energy
break of CR spectrum as well as the corresponding gamma-ray spectrum. As
estimated in Krumholz et al. (2020) the break energy scale as $E_{\rm
br}\sim\left(\frac{\epsilon}{n^{1.5}}\right)^{\frac{1}{\gamma-1}}$, where
$\epsilon$ is the CR energy density, $n$ is the ambient gas density and
$\gamma$ is the spectral index of the CR (assuming a power law distribution).
Thus, $E_{\rm br}$ can be different in different environment. Dramatically, as
we mentioned above, in the case of continuous injection the CR energy density
scale as $\frac{1}{r}$, which predicted a decrease of $E_{\rm br}$ with
increasing $r$. For NGC 6618, such a break are not observed. But for Cygnus
cocoon, a spectral softening at several TeVs have already been detected by
HAWC (Abeysekara et al., 2021), if such a softening is a consequence of the
propagation effects we discussed here, a spatial dependence of the energy
break is expected, which may be tested by later observations.

Table 6:  SED fit results for different radiation models of Src A and Src B

| Model | Index |  $W_{\rm p}$ or $W_{\rm e}$ (erg) a
---|---|---|---
Src A | PP | $3.13^{+0.14}_{-0.12}$ | $4.67^{+0.21}_{-0.18}\times 10^{47}$
IC | $4.35^{+0.18}_{-0.31}$ | $4.63^{+2.31}_{-2.16}\times 10^{50}$
Src B | PP | $\mathrm{2.50_{-0.10}^{+0.10}}$ | $\mathrm{1.83_{-0.10}^{+0.12}\times 10^{45}}$
IC | 3.38${}_{-0.16}^{+0.16}$ | $1.43_{-0.59}^{+0.92}\times 10^{49}$

a The total energy of the parent protons or electrons (with energy $>1$ GeV)
for generate the gamma-ray emission, assuming the medium density
$n(\mathrm{H})$ is $1.3\times 10^{3}\,{\rm cm}^{-3}$ for Src A and $4.4\times
10^{4}\,{\rm cm}^{-3}$ for Src B, and the distance of both Src A and Src B is
${2\,{\rm kpc}}$.

##  5 Conclusions

In this work, we report the detection of gamma-ray emission toward M17 region,
which hosts one of the youngest YMCs in our Galaxy, NGC 6618. In addition to
Cygnus cocoon (Ackermann et al., 2011; Aharonian et al., 2019), NGC 3603 Yang
& Aharonian (2017), W40 (Sun et al., 2020b), and RSGC 1 (Sun et al., 2020a),
this detection provides another case in the population of gamma-ray emitting
YMCs. However, the gamma-ray emission of this source has its unique
characteristic: the gamma-rays can be separated into two components, one is
point-like and with harder spectrum and another is extended and with softer
spectrum. Both components show spatial correlation with dense gases. Assuming
the gamma-rays are produced by the interactions of CRs injected by NGC 6618
with the ambient gas, the different spectral features can be addressed by the
propagation effects investigated in the environments of starburst galaxies
(Krumholz et al., 2020). In this scenario, the MHD turbulent in the ISM can be
damped out effectively, and CRs will stream along field lines and transport
via a process of field line random walk. As a result, the CR diffusion
coefficient is energy-independent at low energy and will become energy-
dependent above certain break energy. The high densities of both gas and CRs
in YMCs are similar to those in starburst galaxies, thus such a comparison is
feasible. In M17 region, the hard spectrum of the point-like source Src B can
be explained by the energy-independent diffusion. As for Src A, due to the
lower gas density and CR density, the damping mechanism can be switched off,
meanwhile, the diffusion coefficient are still the same as in the ISM. The
features related to these propagation effects have not been observed in other
YMC systems yet. But the observed energy break in Cygnus cocoon region can be
explained naturally in such a scenario. Since YMCs can potentially be an
alternative CR source, the propagation of CRs in the vicinity of these sources
would be crucial to understand the injection of CRs to the ISM, which can be
revealed by further observations with more spatial and spectral information.

##  6 Acknowledgements

This research made use of the data from the Milky Way Imaging Scroll Painting
(MWISP) project, which is a multi-line survey in 12CO/13CO/${\rm C}^{18}$O
along the northern galactic plane with PMO-13.7m telescope. MWISP was
sponsored by National Key R&D Program of China with grant 2017YFA0402701 and
CAS Key Research Program of Frontier Sciences with grant QYZDJ-SSW-SLH047.
Bing Liu is supported by the NSFC under grant 12103049. Ruizhi Yang is
supported by the NSFC under grants 12041305, 11421303 and the national youth
thousand talents program in China. Zhiwei Chen acknowledges support from NSFC
grants 11903083, 12173090, U2031202, 11873093, and 11873094.

##  7 Data Availability

The Fermi-LAT data used in this work is publicly available, which is provided
online by the NASA-GSFC Fermi Science Support Center
555https://fermi.gsfc.nasa.gov/ssc/data/access/lat/. A MWISP open data
proposal is required to access the CO line data used in this work, and the
form could be downloaded from http://english.dlh.pmo.cas.cn/op/odp/.

## References

  * Abdo et al. (2009) Abdo A. A., et al., 2009, ApJ, 706, L1
  * Abeysekara et al. (2021) Abeysekara A. U., et al., 2021, Nature Astronomy, 5, 465–471
  * Abramowski et al. (2012) Abramowski A., et al., 2012, A&A, 537, A114
  * Ackermann et al. (2011) Ackermann M., et al., 2011, Science, 334, 1103
  * Aharonian & Atoyan (1996) Aharonian F. A., Atoyan A. M., 1996, A&A, 309, 917
  * Aharonian et al. (2019) Aharonian F., Yang R., de Oña Wilhelmi E., 2019, Nature Astronomy, 3, 561
  * Ballet et al. (2020) Ballet J., Burnett T. H., Digel S. W., Lott B., 2020, arXiv e-prints, p. arXiv:2005.11208
  * Cesarsky & Montmerle (1983) Cesarsky C. J., Montmerle T., 1983, Space Sci. Rev., 36, 173
  * Chen et al. (2012) Chen Z., Jiang Z., Wang Y., Chini R., Tamura M., Nagayama T., Nagata T., Nakajima Y., 2012, PASJ, 64, 110
  * Chen et al. (2013) Chen Z., Nürnberger D. E. A., Chini R., Liu Y., Fang M., Jiang Z., 2013, A&A, 557, A51
  * Chen et al. (2015) Chen Z., Nürnberger D. E. A., Chini R., Jiang Z., Fang M., 2015, A&A, 578, A82
  * Chen et al. (2021) Chen Z., Sun W., Chini R., Haas M., Jiang Z., Chen X., 2021, ApJ, 922, 90
  * Chen et al. (2022) Chen Z., Sefako R., Yang Y., Jiang Z., Yu S., Yin J., 2022, arXiv e-prints, p. arXiv:2204.13296
  * Chibueze et al. (2016) Chibueze J. O., et al., 2016, MNRAS, 460, 1839
  * Chini et al. (1980) Chini R., Elsaesser H., Neckel T., 1980, A&A, 91, 186
  * Domingo-Santamaría & Torres (2006) Domingo-Santamaría E., Torres D. F., 2006, A&A, 448, 613
  * Elmegreen & Lada (1976) Elmegreen B. G., Lada C. J., 1976, AJ, 81, 1089
  * Elmegreen et al. (1979) Elmegreen B. G., Lada C. J., Dickinson D. F., 1979, ApJ, 230, 415
  * Hoffmeister et al. (2008) Hoffmeister V. H., Chini R., Scheyda C. M., Schulze D., Watermann R., Nürnberger D., Vogt N., 2008, ApJ, 686, 310
  * Jiang et al. (2002) Jiang Z., et al., 2002, ApJ, 577, 245
  * Katsuta et al. (2017) Katsuta J., Uchiyama Y., Funk S., 2017, ApJ, 839, 129
  * Krumholz et al. (2020) Krumholz M. R., Crocker R. M., Xu S., Lazarian A., Rosevear M. T., Bedwell-Wilson J., 2020, MNRAS, 493, 2817
  * Lim et al. (2020) Lim W., De Buizer J. M., Radomski J. T., 2020, ApJ, 888, 98
  * Nguyen-Luong et al. (2020) Nguyen-Luong Q., et al., 2020, ApJ, 891, 66
  * Povich & Whitney (2010) Povich M. S., Whitney B. A., 2010, ApJ, 714, L285
  * Povich et al. (2009) Povich M. S., et al., 2009, ApJ, 696, 1278
  * Povich et al. (2016) Povich M. S., Townsley L. K., Robitaille T. P., Broos P. S., Orbin W. T., King R. R., Naylor T., Whitney B. A., 2016, ApJ, 825, 125
  * Ramírez-Tannus et al. (2017) Ramírez-Tannus M. C., et al., 2017, A&A, 604, A78
  * Reid et al. (2019) Reid M. J., et al., 2019, ApJ, 885, 131
  * Saha et al. (2020) Saha L., Domínguez A., Tibaldo L., Marchesi S., Ajello M., Lemoine-Goumard M., López M., 2020, ApJ, 897, 131
  * Stutzki & Guesten (1990) Stutzki J., Guesten R., 1990, ApJ, 356, 513
  * Su et al. (2019) Su Y., et al., 2019, ApJS, 240, 9
  * Sun et al. (2020a) Sun X.-N., Yang R.-Z., Wang X.-Y., 2020a, MNRAS, 494, 3405
  * Sun et al. (2020b) Sun X.-N., Yang R.-Z., Liang Y.-F., Peng F.-K., Zhang H.-M., Wang X.-Y., Aharonian F., 2020b, A&A, 639, A80
  * Wilson et al. (2003) Wilson T. L., Hanson M. M., Muders D., 2003, ApJ, 590, 895
  * Xu et al. (2011) Xu Y., Moscadelli L., Reid M. J., Menten K. M., Zhang B., Zheng X. W., Brunthaler A., 2011, ApJ, 733, 25
  * Yang & Aharonian (2017) Yang R.-z., Aharonian F., 2017, A&A, 600, A107
  * Yang & Wang (2020) Yang R.-Z., Wang Y., 2020, A&A, 640, A60
  * Yang et al. (2018) Yang R.-z., de Oña Wilhelmi E., Aharonian F., 2018, A&A, 611, A77
  * Yin et al. (2022) Yin J., Chen Z., Yao Y., Chen J., Li B., Jiang Z., 2022, Research in Astronomy and Astrophysics, 22, 035021
  * Zabalza (2015) Zabalza V., 2015, in 34th International Cosmic Ray Conference (ICRC2015). p. 922 (arXiv:1509.03319)
  * Zhan et al. (2016) Zhan X.-L., Jiang Z.-B., Chen Z.-W., Zhang M.-M., Song C., 2016, Research in Astronomy and Astrophysics, 16, 56