#  Two eyes, Two views, and finally, One summary! Towards Multi-modal Multi-
tasking Knowledge-Infused Medical Dialogue Summarization

Anisha Saha111Corresponding Author<EMAIL_ADDRESS>Abhisek Tiwari
Sai Ruthvik  Sriparna Saha  Indian Institute of Technology, Patna  Clinical AI
Assistance

###### Abstract

We often summarize a multi-party conversation in two stages: chunking with
homogeneous units and summarizing the chunks. Thus, we hypothesize that there
exists a correlation between homogeneous speaker chunking and overall
summarization tasks. In this work, we investigate the effectiveness of a
multi-faceted approach that simultaneously produces summaries of medical
concerns, doctor impressions, and an overall view. We introduce a multi-modal,
multi-tasking, knowledge-infused medical dialogue summary generation (MMK-
Summation) model, which is incorporated with adapter-based fine-tuning through
a gated mechanism for multi-modal information integration. The model, MMK-
Summation, takes dialogues as input, extracts pertinent external knowledge
based on the context, integrates the knowledge and visual cues from the
dialogues into the textual content, and ultimately generates concise summaries
encompassing medical concerns, doctor impressions, and a comprehensive
overview. The introduced model surpasses multiple baselines and traditional
summarization models across all evaluation metrics (including human
evaluation), which firmly demonstrates the efficacy of the knowledge-guided
multi-tasking, multimodal medical conversation summarization. The code and
dataset are available at https://github.com/NLP-RL/MMK-Summation.

##  1 Overview

Healthcare services are a fundamental necessity for individuals regardless of
their location, social background, or beliefs. Seeking medical guidance is our
typical response to health issues. Nevertheless, the state of healthcare,
especially in government hospitals, is currently distressing in terms of
patient management. In numerous developing nations, waiting several days for a
doctor’s consultation has become the norm due to an inadequate number of
healthcare staff. A recent report from the World Health Organization
anticipates a shortage of 30 million healthcare workers by 2030 relative to
the growing population222https://www.who.int/health-topics/health-workforce.
Moreover, a significant proportion of the population, surpassing one-third,
resides in rural regions333https://www.statista.com/statistics/1328171/rural-
population-worldwide/ where accessing medical professionals remains
challenging. To address these challenges, the utilization of Artificial
Intelligence (AI) based tools and techniques has surged dramatically in the
past few years [2, 4]. One such manifestation is automated symptom
investigation and diagnosis (Figure 1) to assist doctors in conducting
preliminary assessments.

Within a digital healthcare framework, a patient interacts with various
stakeholders, each seeking distinct aspects of the doctor-patient dialogue:
some require the patient’s primary concern, while others seek the doctor’s
perspective or an overarching summary. In this regard, we bring forth a few
concepts, starting with Medical Concern Summary (MCS), which aims to create a
concise overview of the patient’s key issue. Doctor Impression (DI)
encapsulates the concluding reaction and impression of the doctor subsequent
to the conversation with the patient. In the process of clinical diagnosis and
treatment, a patient’s journey typically involves multiple interactions rather
than a single visit. Consequently, reviewing the entire transcript of a
previous lengthy conversation can be time-consuming. Therefore, having access
to the patient’s medical concern summary (MCS) along with the doctor’s
impression (DI) serves as a helpful synopsis/action point of the case for
different healthcare stakeholders, reducing the need to refer to the lengthy
transcript. An example of both MCS and DI are illustrated in Figure 1.
Motivated by the efficacy and importance of MCS and DI, we attempt to study
their impact on overall summary generation and model this information as a
crucial insight while generating an overall summary of discussion between
clinicians and patients.

In our daily lives, the act of summarizing documents is a common task, and the
method employed can vary depending on the context. Nonetheless, a general
approach that many of us tend to follow does exist. As humans, our strategy
for summarizing a dyadic conversation diverges from a strictly sequential
process. We initiate this process by grasping the broader context, after which
we assimilate pertinent specifics from diverse utterances. This enables us to
craft a succinct summary that remains coherent with the overarching context.
Following this initial overview, we take into account the perspectives of each
speaker, which leads us to formulate a comprehensive summary that effectively
captures the central elements. For example, when summarizing the dialogue
(Figure 1), the initial step involves reading the entire conversation.
Subsequently, we identify the central patient concern, referred to as the
medical concern summary (MCS), along with the doctor’s response, termed the
doctor’s impression (DI). Ultimately, the summary of the dialogue is crafted,
focusing on essential elements (MCS and DI), resulting in a coherent and
effective summary. Motivated by the observation, we aim to investigate some
fundamental research questions related to multi-modal medical dialogue
summarization and propose a two-phase, knowledge-infused multi-modal medical
dialogue summarization framework.

Figure 1: Importance of clinical conversation summarization and the two-view
approach for producing an effective summary

Research Questions In this paper, we aim to investigate the following three
research questions: RQ1: Is there any correlation between medical concern
summary (MCS) generation and overall dialogue summarization?, RQ2: Is there
any correlation between doctor impression (DI) generation and overall dialogue
summarization?, RQ3: Are these three tasks correlated? Would learning all
three simultaneously result in an improved overall summary generation?

Table 1: Performances of different baselines and proposed models for multi-
modal clinical conversation summary generation

Model | B-1 | B-2 | B-3 | B-4 | BLEU | R- 1 | R- 2 | ROUGE- L | METEOR | Jaccard Sim | BERT Score
---|---|---|---|---|---|---|---|---|---|---|---
GPT-2 [1] | 11.65 | 5.34 | 2.22 | 0.80 | 5.00 | 21.23 | 4.64 | 20.37 | 23.41 | 0.0717 | 0.6660
BART [3] | 9.94 | 7.18 | 5.16 | 3.72 | 6.50 | 38.37 | 18.04 | 38.50 | 18.68 | 0.1833 | 0.8378
T5 [5] | 42.27 | 32.00 | 24.6 | 18.58 | 29.36 | 54.16 | 31.74 | 51.24 | 43.22 | 0.2582 | 0.8841
MM-MDS | 47.31 | 35.56 | 26.73 | 20.32 | 32.48 | 58.59 | 35.68 | 49.06 | 54.67 | 0.2571 | 0.9081
MM-MDS with MCS | 47.62 | 35.67 | 26.46 | 20.10 | 32.46 | 59.62 | 36.51 | 49.87 | 57.50 | 0.2644 | 0.9115
MM-MDS with DI | 47.05 | 35.27 | 26.42 | 19.79 | 32.13 | 59.46 | 36.29 | 49.76 | 57.11 | 0.2657 | 0.9152
MM-MDS with MCS and DI | 47.70 | 36.48 | 27.52 | 20.85 | 33.14 | 59.74 | 37.59 | 50.93 | 58.85 | 0.2773 | 0.9144
MMK-Summation with MCS | 47.57 | 35.86 | 26.98 | 20.31 | 32.68 | 59.95 | 36.49 | 50.21 | 58.03 | 0.2697 | 0.9145
MMK-Summation with DI | 46.50 | 35.16 | 26.66 | 19.92 | 32.06 | 59.75 | 37.14 | 50.83 | 60.21 | 0.2746 | 0.9162
MMK-Summation | 48.68 | 36.85 | 27.92 | 21.50 | 33.47 | 60.86 | 37.43 | 51.05 | 58.32 | 0.2746 | 0.9180

Table 2: Performances of different baselines and proposed models for medical
concern summary (MCS) generation

Model | B-1 | B-2 | B-3 | B-4 | BLEU | R- 1 | R- 2 | ROUGE- L | METEOR | Jaccard Sim | BERT Score
---|---|---|---|---|---|---|---|---|---|---|---
MM-MDS with only MCS | 29.42 | 14.97 | 1.73 | 0.16 | 11.57 | 46.11 | 20.28 | 41.1 | 30.06 | 0.1958 | 0.6365
MM-MDS | 30.82 | 15.46 | 1.72 | 3.41 | 12.85 | 43.93 | 21.35 | 39.93 | 28.37 | 0.1984 | 0.5811
MMK-Summation with only MCS | 29.53 | 17.02 | 5.35 | 2.92 | 13.71 | 46.49 | 22.87 | 42.88 | 33.69 | 0.1899 | 0.6403
MMK-Summation | 34.77 | 23.92 | 10.79 | 6.90 | 19.10 | 47.86 | 26.67 | 45.18 | 36.83 | 0.2295 | 0.6650

Key Contributions Acquiring insights into medical concerns can facilitate
doctors in formulating appropriate impressions, and the combined knowledge can
contribute to crafting a comprehensive summary. Recognizing the inter-
relationship among these tasks, we build a multi-tasking, multimodal medical
summary generation framework that cohesively learns from all three tasks. With
the immense growth of language models and the proven efficacy of huge pre-
training, we leverage a pre-trained module as the backbone, incorporated with
an adapter-based modality and knowledge infusion. The model receives dialogues
as input, retrieves relevant external knowledge contextualized within the
discourse, combines these knowledge and visual cues expressed in dialogues
with the dialogue text representation, and subsequently produces MCS, DI, and
overall summary. The main contributions of the work are threefold, which are
enumerated below.

  * •

The work first investigates the significance of employing two-stage methods
(homogeneous chunking followed by overall summarization) and multi-tasking
speaker-driven approaches for medical summarization to generate a
comprehensive dialogue summary.

  * •

The proposed MMK-Summation model demonstrates superior performance compared to
numerous baseline models and state-of-the-art approaches, exhibiting a
substantial improvement across various evaluation metrics, including human
assessment.

##  2 Findings

Based on the experiments, we report the following answers (with evidence) and
observations to our investigated research questions (RQ).

RQ1: Is there any correlation between medical concern summary (MCS) generation
and overall dialogue summarization? To address this research question, we
conducted experiments using two models: one that generates an overall summary
solely based on a dialogue (MM-MDS, as shown in Table 1), and another that
simultaneously learns medical summary generation and overall summary
generation tasks (MM-MDS with MCS and MMK-Summation with MCS, also shown in
Table 1). Models that jointly learn these two tasks significantly outperform
the traditional model. The improvements observed across various evaluation
metrics (BELU: 0.20 $\uparrow$, R1-1: 1.36 $\uparrow$, R-2: 0.81 $\uparrow$,
R-L: 1.15 $\uparrow$, METEOR: 3.36 $\uparrow$, Jaccard Sim: 0.013 $\uparrow$
and BERT Score: 0.006 $\uparrow$) with MMK-Summation with MCS provide clear
evidence of the benefits of enhancing overall summary generation through
multitasking with medical concern summary generation.

RQ2: What correlation exists between doctor impression (DI) generation and
overall dialogue summarization? In clinical conversations, two primary
elements stand out: patient utterances and doctor utterances. The precise
summarization of patient utterances aims to capture medical concerns while
summarizing doctor utterances, especially the last few, reveals the doctor’s
impression. Our research revealed a positive correlation between the quality
of medical concern summaries and overall summarization, and vice versa.
Consequently, we inferred that the doctor’s impression also influences the
overall summary. Interestingly, our findings indicate that the doctor’s
impression alone doesn’t significantly contribute to the overall summary
compared to using medical concern summaries alone (MM-MDS versus MMK-Summation
with MCS). However, when combined with medical concern summaries, it notably
enhances the performance of overall summary generation (MMK-Summation with MCS
versus MMK-Summation). Upon delving into the rationale, we observed that
medical concern summaries exhibit a stronger correlation with the overall
summary compared to doctor impressions, evident in both word overlap and
semantic similarity. Thus, based on this evidence, we assert a minor yet
discernible impact of doctor impressions on overall summary generation.

RQ3: Are these three tasks correlated? Would learning all three simultaneously
result in an improved overall summary generation? A dialogue summary is a
cohesive piece of concise information encapsulating key points from different
speakers engaged in the discussion. Thus, we, as humans, summarize a
conversation by first summarizing each speaker’s context and then forming a
unified dialogue summary. Humans do it very fast, and thus, it seems like it
is a one-step process. Motivated by this, we build a multitasking model, MMK-
Summation, which takes dialogue as input and generates MCS, DI, and overall
summary simultaneously. The proposed multitasking MMK-Summation model
outperforms both MM-MDS (generating only the overall summary) and MMK-
Summation with MCS/DI. The obtained improvement across different metrics are
as follows: BLEU: $0.99\uparrow$, R-1: $2.27\uparrow$, R-2: $1.75\uparrow$,
ROUGE-L: $1.99\uparrow$, METEOR: $3.65\uparrow$, Jaccard Sim: $0.017\uparrow$,
and BERT Score: $0.009\uparrow$. The improvements in terms of BLEU, ROUGE,
METEOR, and BERT Score are significant for a generation task, establishing
confidence in our hypothesis that learning these three tasks, MCS, DI, and
overall summary generation tasks, would help in yielding an efficient overall
summary.

Key Observations The key observations and findings are as follows: (i) The
MMK-Summation model, which produces MCS, DI, and overall summaries from the
decoder, surpasses the MM-MDS model with MCS/DI, which generates MCS & DI from
the encoder and overall summary from the decoder. This discovery suggests that
context-aware cross-attention in the decoder notably impacts (improves)
context embedding necessary for effective decoding. Additionally, it aligns
with human behavior, indicating instantaneous speaker-driven chunking followed
by summarization. (ii) The task of generating medical concern summaries is
more closely correlated with overall summary generation than with doctor
impression. (iii) We noted that conventional models like BART and T5 are
notably deficient in two main areas: (a) comprehending visuals and (b)
maintaining symptom-consistent diagnoses. Thus, incorporating contextualized
M-modality fusion of visuals and knowledge has led to significantly superior
performance in terms of all evaluation metrics, including human evaluation.

## References

  * Budzianowski and Vulic [2019] P. Budzianowski and I. Vulic.  Hello, it’s gpt-2-how can i help you? towards the use of pretrained language models for task-oriented dialogue systems.  _EMNLP-IJCNLP 2019_ , page 15, 2019.
  * Chang et al. [2019] Y.-W. Chang, P.-Y. Hsu, Y. Wang, and P.-Y. Chang.  Integration of online and offline health services: The role of doctor-patient online interaction.  _Patient Education and Counseling_ , 102(10):1905–1910, 2019.
  * Lewis et al. [2019] M. Lewis, Y. Liu, N. Goyal, M. Ghazvininejad, A. Mohamed, O. Levy, V. Stoyanov, and L. Zettlemoyer.  Bart: Denoising sequence-to-sequence pre-training for natural language generation, translation, and comprehension.  _arXiv preprint arXiv:1910.13461_ , 2019.
  * Lin et al. [2021] S. Lin, P. Zhou, X. Liang, J. Tang, R. Zhao, Z. Chen, and L. Lin.  Graph-evolving meta-learning for low-resource medical dialogue generation.  In _Proceedings of the 35th AAAI Conference on Artificial Intelligence_ , pages 13362–13370. AAAI Press, 2021.
  * Raffel et al. [2020] C. Raffel, N. Shazeer, A. Roberts, K. Lee, S. Narang, M. Matena, Y. Zhou, W. Li, and P. J. Liu.  Exploring the limits of transfer learning with a unified text-to-text transformer.  _The Journal of Machine Learning Research_ , 21(1):5485–5551, 2020.


#  Proactive and Resilient UAV Orchestration for QoS Driven Connectivity and
Coverage of Ground Users

Yuhui Wang and Junaid Farooq
Department of Electrical and Computer Engineering,
University of Michigan-Dearborn, Dearborn, MI 48128 USA,
Emails: {ywangdq<EMAIL_ADDRESS>
###### Abstract

Unmanned aerial vehicles (UAVs) are being successfully used to deliver
communication services in applications such as extending the coverage of 5G
cellular networks in remote areas, emergency situations, and enhancing the
service quality in regions of dense user populations. While optimized
placement solutions have been proposed in literature for ensuring quality of
service (QoS) of users, they may not be ideal in highly mobile, autonomous,
and diverse network scenarios. This paper proposes a proactive and resilient
framework for distributed and dynamic orchestration of UAV small cells to
provide QoS differentiation in the network. The UAV locations are tailored to
end-user locations and service needs while ensuring that UAVs maintain
localized backhaul connectivity. Simulation experiments show that under
scenarios where physical placement can achieve service differentiation, the
developed framework leads to a stable configurations of UAVs satisfying above
90% of user QoS requirements.

###### Index Terms:

Unmanned Aerial Vehicles, Connectivity, Resilience, Quality-of-Service,
Swarming.

##  I Introduction

UAVs have proven to be a convenient choice for providing on-demand wireless
connectivity in disasters and emergency scenarios [1]. They are also critical
in enhancing and improving the service gaps in fifth generation (5G) cellular
networks. In certain situations such as public events, gatherings, parades, or
sporting events, there may be a need to supplement the cellular networks with
small cell base stations (BSs) mounted on UAVs [2, 3]. While several factors
play a role in the quality of communication experienced by the users connected
to the UAVs, the placement the UAV BSs can be crucial in ensuring the strength
of wireless signals between the UAV and the ground users [4]. Moreover,
multiple types of users may be present in the network, e.g., those that have
subscribed for higher speeds and preferred resource allocation, hereby
referred to as _premium users_ , and those that are using standard service,
hereby referred to as _regular users_. Therefore, the placement of UAVs needs
to be optimized to cater for the quality of service (QoS) requirements of
different users [5, 6]. An illustration of the considered network scenario is
shown in Fig. 1, where the premium user and regular user are connected to a
network of UAV BSs that are orchestrated in the air to meet the data needs of
the users.

In conventional cellular networks, the BS locations are fixed and there is not
much that can be done to improve the link quality of a particular user from a
path-loss or interference standpoint. However, in the case of a mobile aerial
BS, there is a high degree of freedom in the location of the BS, allowing it
to move closer to a particular user and vice versa to satisfy QoS requirements
[7]. For instance, if the user has a premium subscription and if it does not
achieve sufficient data rates due to interference or other channel
impairments, the network can reconfigure itself to improve the QoS of the
user. In essence, we can leverage the mobility and flexibility of aerial
placement in UAVs to provide a customized service to ground users. This
enables new service models and pricing options by the network operators that
will trigger higher revenue generation. Existing cellular operators can also
extend their network by adding on-demand small cells in scenarios where a
higher user density is expected [8].

Figure 1: Orchestration of UAVs in the air to serve premium and regular ground
users with different QoS requirements.

There are several existing works in the literature that consider the QoS
provisioning in the placement and resource allocation decisions of UAVs [9,
10, 11, 12, 13, 14, 15]. These consider factors such as energy efficiency [14,
16], QoS guarantees [15, 9, 10], and trajectory optimization [11, 17, 12].
However, there are several limitations: (i) The existing solutions are
centralized and are not based on feedback from the users, whereas our proposed
framework proactively involves the performance of users into the placement
decisions via real-time feedback on the achieved data rates. (ii) They are not
adaptive to changes in user locations and performance needs over time, whereas
our proposed solution is dynamic and continuously adapts to the highly mobile,
autonomous, and dynamic network scenarios. Hence it is resilient to potential
failures and cyber-physical attacks on the network. (iii) Existing works focus
on the placement decisions while ignoring the details of UAV dynamics and
control aspects, whereas we have used a dynamical model to realistically mimic
the practical scenarios and used a controller based approach inspired from
swarming and flocking of robotic systems.

Although the problem of optimized placement of UAVs can be solved in a
centralized manner, it becomes less applicable in complex real world scenarios
where a centralized planner may not be available to have a global view of the
network. We propose a completely distributed and dynamic approach to
orchestrate the UAVs using an incentive mechanism design approach that
facilitates the fulfilment of QoS requirements of the users. Moreover, our
framework can also be augmented to centralized placement solutions to further
fine tune the performance of users in the network and adapt to their changing
positions and QoS needs. The key contributions of this work are as follows:

  * •

We develop a proactive and resilient dynamic framework for real-time UAV BS
orchestration that provides QoS differentiation for premium and regular ground
users.

  * •

We propose a UAV swarming approach that ensures localized backhaul
connectivity among UAVs while adaptively orienting themselves to cater for
user QoS requirements.

  * •

We conduct numerical experiments in several simulation scenarios to test our
solution. The results demonstrate above 90% efficacy of the proposed
algorithm.

The rest of the paper is organized as follows: Section II provides an overview
of the system model and defines the QoS provisioning problem, Section III
details the methodology of our proposed approach, Section IV provides
simulation based validation of our framework, followed by the conclusion in
Section V.

##  II System Model

Consider a set of ground users $\mathcal{M}=\\{1,\ldots,M\\}$ and a set of
UAVs $\mathcal{U}=\\{1,\ldots,U\\}$, equipped with small cell BSs, to provide
wireless connectivity to the ground users. The Cartesian coordinates of ground
users at time $t$ are denoted by
$\boldsymbol{y}(t)=[y_{1}(t),y_{2}(t),\ldots,y_{M}(t)]^{T}$, where
$y_{i}(t)\in\mathbb{R}^{3},\forall i\in\mathcal{M},t\geq 0$. Note that the
height of all ground users is assumed to be zero. Similarly, the Cartesian
coordinates of the UAVs at time $t$ are denoted by
$\boldsymbol{q}(t)=[q_{1}(t),q_{2}(t),...,q_{U}(t)]^{T}$, where
$q_{i}(t)\in\mathbb{R}^{3},\forall i\in\mathcal{U},t\geq 0$. The height of all
the UAVs is assumed to be $H$111The altitudes of UAVs is generally regulated
based on the flying environment and safety standards. Therefore, we assume a
predefined constant elevation for all UAVs.. The velocity of the UAVs at time
$t$ is denoted by $\boldsymbol{p}(t)=[p_{1}(t),p_{2}(t),...,p_{U}(t)]^{T}$,
where $p_{i}(t)\in\mathbb{R}^{3},\forall i\in\mathcal{U},t\geq 0$. Since the
UAVs maintain a fixed height $H$, the vertical velocity component of all UAVs
is zero. The UAVs have a maximum communication range of $r\in\mathbb{R}^{+}$
such that any two UAVs can communicate only if the Euclidean distance between
them is less than $r$. The neighboring UAVs of a typical UAV denoted by $i$ is
represented by the set $N_{i}=\\{j\in\mathcal{U}\backslash
i:\|q_{i}-q_{j}\|\leq r\\},\forall i\in\mathcal{U}$. User $m$ connects to UAV
$i$ if $i=\underset{k}{\arg\min}\|q_{k}-y_{m}\|,\forall k\in\mathcal{U}$. The
set of users served by UAV $i$ is represented by $\mathcal{M}_{i}^{C}$, which
is limited by the maximum serving capacity $|\mathcal{M}_{i}^{C}|\leq
N^{\max}$.

Each UAV uses a wireless channel $\mathcal{L}_{i},i\in\mathcal{U}$ and
$\mathcal{L}_{i}\in\\{0,1,2,3,...\\}$ to communicate with ground users. We
consider two types of users in the network denoted by $\chi$, where
$\chi(m)=\text{(prem)}$ if user $m$ is a premium user and
$\chi(m)=\text{(reg)}$ if user $m$ is a regular user. Regular users require a
data rate $\rho^{(\text{reg})}$ while premium users require a data rate
$\rho^{(\text{prem})}$. Note that $\rho^{(\text{prem})}>\rho^{(\text{reg})}$
signifying a higher data rate requirement for premium users. Regular users can
only connect to UAVs through the default channel $L_{0}$, while premium users
can access all channels $L_{i}$, $i\geq 0$. In the default configuration, all
UAVs use channel $L_{0}$ to communicate with users, i.e.,
$\mathcal{L}_{i}=L_{0},\forall i\in\mathcal{U}$. If a premium user cannot
achieve the target data rates, the connected UAV may switch channel to reduce
interference from other UAVs and consequently improve performance.

###  II-A Path-Loss Model

The wireless signal from the UAV to the ground users may encounter a line-of-
sight (LoS) path or a non-line-of-sight (NLoS) path due to reflections from
buildings and other obstacles. The signal attenuation in each case is
different and we use a probabilistic approach to model the availability of
these paths. The probability of having a LoS path between the user and UAV is
denoted by $p_{\text{LoS}}$ and the probability of having an NLoS path is
denoted by $p_{\text{NLoS}}$ with $p_{\text{NLoS}}=1-p_{\text{LoS}}$. The
probability of LoS connection is computed as [18]:

$p_{\text{LoS}}=\frac{1}{1+\vartheta\exp\left(-\xi\frac{180}{\pi}\phi-\vartheta\right)},$
(1)

where $\vartheta$ and $\xi$ are constants depending on the environment and
$\phi$ is the elevation angle as shown in Fig. 2.

Figure 2: Relative positions of UAV and ground user.

NLoS paths may experience higher path loss as compared to LoS connections. We
use the widely accepted air-to-ground path-loss model in [18] to characterize
the path-loss between UAVs and ground users. The path loss between a UAV $i$
and user $m$ can be expressed as:

$\displaystyle PL_{i,m}=$ $\displaystyle 10\log_{10}\left(\frac{4\pi
f_{c}D_{i,m}}{c}\right)^{\delta}+p_{\text{LoS}}\hskip
2.168pt\eta_{\text{LoS}}+p_{\text{NLoS}}\hskip 2.168pt\eta_{\text{NLoS}},$ (2)

where $D_{i,m}$ is the distance between the UAV and ground user, i.e.,
$D_{i,m}=\|q_{i}-y_{m}\|$, $i\in\mathcal{U},m\in\mathcal{M}$, $c$ denotes the
speed of light, $f_{c}$ is the carrier frequency, and $\delta$ is the path-
loss exponent. The first term in (2) represents the free-space path-loss,
while $\eta_{\text{LoS}}$ and $\eta_{\text{NLoS}}$ are the additional average
losses for LoS and NLoS paths respectively.

###  II-B Quality of Service (QoS)

To quantify the QoS of ground users, we use the Shannon capacity based on the
signal-to-interference-plus-noise ratio (SINR). The downlink SINR of UAV $i$
connecting to user $m$ can be expressed as:

$\displaystyle\text{SINR}_{i,m}$
$\displaystyle=\frac{P_{r}^{i,m}}{\sigma^{2}+{\sum_{k\in\mathcal{U}\backslash
i,\mathcal{L}_{i}=\mathcal{L}_{k}}P_{r}^{k,m}}},$ (3)

with $P_{r}^{i,m}=P_{t}^{i,m}/10^{PL_{i,m}/10}$. $P_{r}^{i,m}$ is the received
signal power by the user and $P_{t}^{i,m}$ is the transmission power of the
UAVs. $\sum_{k\in\mathcal{U}\backslash
i,\mathcal{L}_{i}=\mathcal{L}_{k}}P_{r}^{k,m}$ is interference power from
other UAVs using the same channel $\mathcal{L}_{i}=\mathcal{L}_{k}$, and
$\sigma^{2}$ is the thermal noise power. $P_{r}^{i,m},P_{t}^{i,m},P_{r}^{k,m}$
and $\sigma^{2}$ are all in the units of mW. Consequently, the achieved data
rate of UAV $i$ connecting to user $m$, denoted by $\mathcal{C}_{i,m}$ can be
written as follows:

$\mathcal{C}_{i,m}=B\hskip 1.4457pt\log_{2}(1+\text{SINR}_{i,m}),$ (4)

where $B$ is the available channel bandwidth. Since user $m$ is always
connected by its closest UAV, we subsequently drop the subscript $i$ for
brevity of notation and use $\mathcal{C}_{m}$ to denote the QoS achieved by
user $m$ while served by its associated UAV.

###  II-C Problem Formulation

The objective of the proposed framework is to determine a set of positions at
each time instant that minimizes the difference between achieved QoS,
$\mathcal{C}_{m}$ and target QoS, $\rho^{\chi(m)}$ for all users
$m\in\mathcal{M}$, while ensuring that the spatial interval of UAVs to be
between the minimum safe distance $d$ and the maximum communication range $r$.
This optimization problem can be formulated as follows:

$\mathbf{P0}:\hskip
28.90755pt\min_{\textbf{q}(t)}\sum_{m\in\mathcal{M}}|\mathcal{C}_{m}-\rho^{\chi(m)}|,$
(5) s.t. $\displaystyle\|q_{i}(t)-y_{j}(t)\|\leq r,\hskip 7.22743pt\forall
i\in\mathcal{U},j\in\mathcal{M}_{i}^{C},$ (6) $\displaystyle
d\leq\|q_{i}(t)-q_{j}(t)\|\leq r,\hskip 7.22743pt\forall i\in\mathcal{U},j\in
N_{i}.$

In problem $\mathbf{P0}$, we have a set of facilities (UAV $\mathcal{U}$) and
a set of clients (users $\mathcal{M}$). Each user has a demand of target data
rate $\rho^{\chi(m)}$ and is served by one UAV. The constraints in (6) signify
that the Euclidean distance between any pair of UAV and its connected user,
i.e., $\|q_{i}(t)-y_{j}(t)\|$, should be less than the maximum communication
range $r$ of the UAVs. In addition, the Euclidean distance between UAVs, i.e.,
$\|q_{i}(t)-q_{j}(t)\|$, should also in the the range of $(d,r)$ to avoid
collision and ensure connectivity. This problem is a special case of the
classical capacitated facility location problem (CFLP) which is known to be
NP-hard [19]. Therefore, the problem P0 is also an NP-hard problem and we need
to resort to heuristic approaches to obtain acceptable solutions.

##  III Methodology

In this section, we use a dynamic and distributed approach to solve the
optimization problem described by P0. The UAV network is orchestrated in real-
time according to the QoS needs of the users through a controller-based UAV
swarming approach, which autonomously creates a resilient configuration and
leverages channel switching to provide customized service to both regular and
premium users. We use the widely accepted kinematic model in robotics and
control literature to describe the dynamics of the UAVs as
$\dot{q_{i}}=p_{i}$, and $\dot{p_{i}}=u_{i}$, where $u_{i}$ is the control
input, $q_{i},p_{i},u_{i}\in\mathbb{R}^{3}$ and $i\in\mathcal{U}$. For the
UAVs to organize themselves in a distributed manner, the control input $u_{i}$
is designed as a linear combination of the following three terms:

$u_{i}=f_{i}+g_{i}+h_{i},$ (7)

where $f_{i}$ is an inter-UAV attractive/repulsive term, $g_{i}$ is a velocity
consensus term, $h_{i}$ is a term defining the individual goal of each UAV. In
order to create smooth incentive functions with finite cut-offs, we leverage
the bump function $\alpha_{\\{\gamma\\}}(z)$ in [20] and $\sigma$-norm in
[21]. The bump function is defined as follows:

$\alpha_{\\{\gamma\\}}(z)=\left\\{\begin{array}[]{ll}1,&\text{if}\ 0\leq
z<\gamma,\\\
\frac{1}{2}\left(1+\cos(\pi\frac{z-\gamma}{1-\gamma})\right),&\text{if}\
\gamma\leq z<1,\\\ 0,&\text{if}\ z\geq 1,\end{array}\right.$ (8)

where $\gamma\in[0,1]$ and the $\sigma$-norm is defined as
$\|z\|_{\sigma}=\frac{1}{\epsilon}\left(\sqrt{1+\epsilon\|z\|^{2}}-1\right)$,
where $\epsilon>0$ is a positive constant. The gradient of the $\sigma$-norm
can be expressed as
$\nabla\|z\|_{\sigma}=\frac{z}{\sqrt{1+\epsilon\|z\|^{2}}}=\frac{z}{1+\epsilon\|z\|_{\sigma}}$.
We make use of these function to create incentives for the UAVs to achieve the
optimal solution to P0. We describe each of these incentive functions along
with their role in the following subsections.

###  III-A UAV Swarming

The swarming behavior requires UAVs to maintain a certain distance with each
other and follow a group objective. This is achieved by introducing the
attraction/repulsion and velocity consensus functions, which maintains the
distance between UAVs through the control of their velocity. Considering the
objective to enhance QoS for users, a UAV $i$ should have the tendency to
approach another UAV $j$ if UAV $j$ exceeds its serving capacity. Thus, the
attractive/repulsive function can be written as follows [22, 23]:

$f_{i}(\mathbf{q},\mathcal{M}^{c})=\sum_{j\in
N_{i}}\Bigg{[}\Phi(\|q_{j}-q_{i}\|_{\sigma})+\\\
a\left(1-\alpha_{\\{0\\}}\left(\frac{\|(|\mathcal{M}^{c}_{i}|-N^{\max})^{+}\|_{\sigma}}{\|N^{\max}\|_{\sigma}}\right)\right)\Bigg{]}\textbf{v}_{i,j},$
(9)

where $\textbf{v}_{i,j}=\nabla\|q_{j}-q_{i}\|_{\sigma}$ is the vector from
$q_{i}$ to $q_{j}$. The function $\Phi(z)$ is expressed as:

$\Phi(z)=\alpha_{\\{0.2\\}}\left(\frac{z}{\|r\|_{\sigma}}\right)\phi(z-\|d\|_{\sigma}),$
(10)

where $\phi(z)=\frac{1}{2}[(a+b)\frac{(z+c)}{\sqrt{1+(z+c)^{2}}}+(a-b)]$ and
$c=|a-b|/\sqrt{4ab}$ to ensure that $\phi(0)=0$. Here, $r$ is the maximum
communication range and $d$ is the minimum distance between UAVs. The velocity
consensus function works as a damping force that leads to a match of movements
between neighboring UAVs. The function $g_{i}(\mathbf{p},\mathcal{M}^{c})$ is
in the form:

$g_{i}(\mathbf{p},\mathcal{M}^{c})=\sum_{j\in
N_{i}}\alpha_{\\{0.2\\}}\left(\frac{\|q_{i}-q_{j}\|_{\sigma}}{\|r\|_{\sigma}}\right)(p_{j}-p_{i}).$
(11)

This enforces the alignment between UAVs to increase as the distance between
them decreases, until reaching the maximum value at a distance of $0.2r$. The
velocity consensus function reduces potential collisions and disconnections
between UAVs during real-time operations.

Algorithm 1 QoS-Driven Dynamic Orchestration

1:Initialize position, velocity, and channel for each UAV $q_{i}(0)$,
$p_{i}(0)$, $\mathcal{L}_{i}(0)\leftarrow L_{0}$; Assign user category
$\chi(m)$ and target data rate $\rho^{\chi(m)}$, record mean data rates
$\mu_{m}$ in last $\tau$ seconds for each user $m\in\mathcal{M}$.

2:while target data rates not satisfied do

3: for each UAV $i$ do

4: Determine the set of connected users $\mathcal{M}_{i}^{C}$.

5: for user $m\in\mathcal{M}_{i}^{C}$ do

6: Compute the SINR and data rates $\mathcal{C}_{m}$ for each

7: user $m$ using (3) and (4).

8: if $\chi=(\text{prem})$ & $\mathcal{C}_{m}<\rho^{\chi(m)}$ &
$\mathcal{C}_{m}\leq\mu_{m}$ then

9: Premium user $m$ and UAV $i$ select a new

10: available channel for communication.

11: end if

12: Recompute mean data rates $\mu_{m}$.

13: end for

14: Compute control input $u_{i}$ for UAV $i$ using (7).

15: Update the position $q_{i}$ and velocity $p_{i}$ using the

16: discretized UAV dynamics.

17: end for

18:end while

###  III-B QoS-Driven Dynamic Orchestration

The QoS requirements of users may vary by their subscribed services and
application scenarios. We define a QoS-driven incentive function which aims to
improve the received data rates of users. Based on target data rates
$\rho^{\chi(m)}$ and achieved data rates $\mathcal{C}_{m}$ for each user
$m\in\mathcal{M}_{i}^{c}$ connected to UAV $i$, this term reflects the
affinity of a UAV and its connected users:

$\displaystyle
h_{i}(\mathbf{q},\mathcal{M}^{c})=\sum_{m\in\mathcal{M}\backslash\mathcal{M}_{i}^{c}}c_{1}\left(\frac{\rho^{\chi(m)}-\mathcal{C}_{m}}{\rho^{\chi(m)}}\right)^{+}\nabla\|q_{i}-q_{m}\|_{\sigma}$
(12)
$\displaystyle+\sum_{m\in\mathcal{M}_{i}^{c}}c_{2}^{(\chi)}\alpha_{\\{0\\}}\left(\frac{\mathcal{C}_{m}}{\beta\rho^{\chi(m)}}\right)\phi(\rho^{\chi(m)}-\mathcal{C}_{m})\nabla\|q_{m}-q_{i}\|_{\sigma}.$

The first part of the equation represents a repulsive force from users against
other UAVs. If a user’s data rate requirement is not satisfied, it may
generate a repulsive force against UAVs that it is not connected to for
reducing interference. The second part of the equation represents an
attractive force from users that are connected with UAV $i$ towards UAV $i$.
If a user’s data rate requirement is not satisfied, it may generate an
attractive force to the UAV it connects with to improve link quality. Note
that the parameter $c_{2}$ is selected differently for regular and premium
users, i.e., $c_{2}^{(\text{prem})}=1.5c_{2}^{(\text{reg})}$.

The framework for QoS-driven dynamic placement of UAVs can be summarized using
Algorithm 1. We assume each user always connect with the nearest UAV from its
current position. In every iteration, UAVs compute data rates for all
connected users and keep records of them. If the target data rates for a
premium user is not reached and the real-time data rates have not increased
during the past $\tau$ seconds, the UAV will switch to a free channel to
reduce interference from other UAVs and therefore increase the data rates. The
control input $u$ is computed using (7) and is used to update position of the
UAVs. The proposed distributed algorithm has linear time complexity of
$\mathcal{O}(\|N\|+\|\mathcal{M}\|)$ for each UAV, where $N$ is the set of
neighboring UAVs and $\mathcal{M}$ is the set of all users within the
communication range of the UAV.

##  IV Simulation Results

In this section, we simulate different UAV network scenarios and validate the
efficacy of our proposed algorithm. Unless otherwise specified, the model
parameters used in the simulations are provided as follows: the height of all
UAVs, $H=100$ m ; the minimum distance between UAVs, $d=100$ m; the maximum
communication range of UAVs, $r=300$ m. The target data rates for premium and
regular users are set to $\rho^{(\text{prem})}=300$ Mbps and
$\rho^{(\text{reg})}=100$ Mbps respectively. Other parameters and their values
are tabulated in Table 1.

Parameter | Notation | Value
---|---|---
UAV carrier frequency | $f_{c}$ | 2 GHz
LoS additional avg loss | $\eta_{LoS}$ | 0.1 dB
NLoS additional avg loss | $\eta_{NLoS}$ | 21.0 dB
Environment parameters | $\vartheta,\xi$ | 4.88, 0.43
UAV transmission power | $P_{t}$ | 37.0 dBm
Bandwidth | $B$ | 15 MHz
Thermal Noise | $\sigma^{2}$ |  $-$80 dBm
Incentive function | $\epsilon,a,b,N^{\max}$ | 0.1, 5.0, 5.0, 80.0
parameters | $c_{1},c_{2}^{(\text{reg})},\beta$ | 6.0, 4.0, 1.5
Channel switch interval | $\tau$ | 5.0 sec
TABLE I: Simulation parameters Figure 3: Simulation results with 3 users (2
premium and 1 regular) and two UAVs using the proposed method. Red lines show
the trajectories of UAVs. All three users achieve close to target data rates
on convergence. Figure 4: Comparison simulation using flocking algorithm only.
The users cannot achieve target data rates.

###  IV-A QoS Driven Proactive UAV Placement

We first consider a simple scenario of three users served by two UAVs to
demonstrate the operation of the proposed framework. The user locations and
types are shown in Fig. 3. We consider two premium users and one regular user.
UAVs are initially placed arbitrarily and the red dotted lines shows how their
trajectory evolve over time to reach the final configuration. At the initial
state, the premium users do not achieve the required data rates so they
generate an attractive force on their connected UAVs. Hence, we observe that
the two UAVs move closer to the respective premium users. However, to further
enhance data rates, the UAVs need to move farther from each other to reduce
interference. Therefore, one of the UAVs follows an arc trajectory around the
connected user to maximize the data rates of the regular user while
maintaining the already achieved data rates of the premium user. Note that the
UAVs also need to ensure a minimum distance between then to enable localized
backhaul connectivity. This is highlighted by the black dotted line in the
figure. Simulation using standard flocking algorithm [21] alone was compared
in Fig. 4. The UAVs can gather around the centroid of the three users and
maintain the minimum distance with each other. However, they cannot adapt
their positions according to the QoS needs and resulted in unsatisfied data
rates.

###  IV-B Resilience under Failures and Attacks

In Fig. 5, we consider the scenario of 600 users located in a rectangular
region mimicking a parade or a procession. The users are placed randomly
according to a uniform distribution and are split into two types. We assume
that 30% of the users are premium users and are located in the left part
highlighted by green dots. The regular users are located on the right side and
highlighted in blue. 15 UAVs are deployed to provide service to the users. The
initial positions and traces of the UAVs are shown by dotted red lines in Fig.
5. Due to the signal interference and geographical distance, the premium users
gain low data rates in the initial state. At time $t=9.0$s, the data rates of
premium users reach a bottleneck so UAVs which serve premium users switch
their channels because the premium users cannot meet their QoS requirements.
Seven UAVs change to serve premium users only and the other three UAVs serve
the regular users. As the interference is eliminated, the mean data rates for
premium users goes up to the desired 300 Mbps at time $t=10.0$s. To test the
resilience of the network, 30% of UAVs randomly fail at t=15s. Failed UAVs are
labeled with red cross. The UAVs rapidly reconfigure their formation and
recover the QoS. At t=30s, the mean data rates for premium users recover to
over 293 Mbps and the data rates for regular users recover to over 79 Mbps.

Figure 5: Simulation results with 600 users randomly distributed in a
rectangular space. 30% of users are premium users and the rest are regular
users. 30% of UAVs which randomly failed at t=15s were labeled with red cross.

###  IV-C Network Size and Impact on QoS Provisioning

Next, we investigate the impact of changing the number of UAV BSs on the the
achieved performance of the ground users. We use the scenario illustrated in
Fig. 5 and increase the number of UAVs to observe the number of served users
and achieved data rates as compared to the targets of both premium and regular
users. Fig. 6(a), 6(b), and 6(c) respectively show the percentage of served
users, average data rates, and percentage of fulfilled data requirements for
different user groups against the number of UAVs. When only 6 UAVs are
available, the UAVs can only serve a small proportion of the total users due
to capacity constraints and can be far away from users that are located at the
edges. So the percentage of served users and the average data rates are low.
As the number of UAVs increases to 21, the QoS of the network improves until
all the premium users are served and achieve the required data rates. It can
be shown that the QoS of the network makes fast improvements at initial states
before converging to their maximum values as the number of UAVs increases
until reaching the capacity limit of the UAV network. Note that the premium
users achieve higher performance not only in terms of the target data rates
but also in the percentage of served users and the percentage of fulfilled
requirements. The only exception to this is observed in Fig. 6(c) for small
number of UAVs since the target data rates of premium users are too aggressive
and could not be achieved by a small number of UAVs. However, with large
number of UAVs, both the primary users and regular users are able to
completely achieve their targets.

(a)

(b)

(c)

Figure 6: Percentage of served users, average data rates and percentage of
fulfilled data rate requirements for different user groups against the number
of UAVs.

##  V Conclusions

In this paper, we investigated the design of QoS-centric dynamic placement of
UAVs to provide customized service quality for different user requirements. We
leveraged the mobility and flexibility of aerial placement in UAVs and
developed on-demand UAV small cell deployment strategies to not only provide
higher data rates to users in required regions but also provide service
differentiation for meeting users needs (live streaming, social media
activities) and also higher revenues for service providers. The proposed
method utilized distributed UAV flocking and swarming dynamics, and extended
it by adding an incentive function based on target and real-time data rates.
Compared with open-loop optimization approaches, the QoS-driven dynamic
placement algorithm has linear time complexity and provides adaptive coverage
and differentiated service levels. The simulation results show that our
proposed method is able to fulfill the data rates requirements of premium and
regular users in both simple and complex scenarios. The proposed solution can
be useful in real networks where fixed small cell deployment may be
economically in-viable for occasional needs for high speed connectivity.

## References

  * [1] Y. Wang and J. Farooq, “Resilient UAV formation for coverage and connectivity of spatially dispersed users,” in _IEEE International Conference on Communications (ICC 2022)_ , Seoul, South Korea, May 2022.
  * [2] “Integration of civil unmanned aircraft systems (UAS) in the national airspace system (NAS) roadmap,” Federal Aviation Administration, Tech. Rep., 3rd Edition, Rep. no. AUS-200807-002, 2020.
  * [3] K. G. Panda, S. Das, and D. Sen, “Efficient UAV placement strategy for guaranteed QoS demand,” in _IEEE 92nd Vehicular Technology Conference (VTC-Fall 2020)_ , Victoria, BC, Canada, Dec. 2020.
  * [4] Y. Zeng, Q. Wu, and R. Zhang, “Accessing from the sky: A tutorial on UAV communications for 5G and beyond,” _Proceedings of the IEEE_ , vol. 107, no. 12, 2019.
  * [5] C. Zhang, L. Zhang, L. Zhu, T. Zhang, Z. Xiao, and X.-G. Xia, “3D deployment of multiple UAV-mounted base stations for UAV communications,” _IEEE Transactions on Communications_ , vol. 69, no. 4, 2021.
  * [6] S. Zhang and W. Cheng, “Statistical QoS provisioning for UAV-enabled emergency communication networks,” in _IEEE Globecom Workshops (GC Wkshps 2019)_ , Big Island, HI, USA, Dec. 2019.
  * [7] J. Luo, J. Song, F. Zheng, L. Gao, and T. Wang, “User-centric UAV deployment and content placement in cache-enabled multi-UAV networks,” _IEEE Transactions on Vehicular Technology_ , pp. 1–1, 2022.
  * [8] H. R. Chi and A. Radwan, “An overview of on-demand deployment optimization of small cells,” _IEEE Network_ , vol. 35, no. 2, 2021.
  * [9] H. Niu, X. Zhao, and J. Li, “3D location and resource allocation optimization for UAV-enabled emergency networks under statistical QoS constraint,” _IEEE Access_ , vol. 9, 2021.
  * [10] X. Zhang, J. Wang, and H. V. Poor, “Statistical delay and error-rate bounded QoS provisioning for 6G mURLLC over AoI-driven and UAV-enabled wireless networks,” in _IEEE Conference on Computer Communications (INFOCOM 2021)_ , 2021, pp. 1–10.
  * [11] Y. Huang, M. Cui, G. Zhang, and W. Chen, “Bandwidth, power and trajectory optimization for UAV base station networks with backhaul and user QoS constraints,” _IEEE Access_ , vol. 8, 2020.
  * [12] L. Shi and S. Xu, “UAV path planning with QoS constraint in device-to-device 5G networks using particle swarm optimization,” _IEEE Access_ , vol. 8, 2020.
  * [13] J. Qin, Z. Wei, C. Qiu, and Z. Feng, “Edge-prior placement algorithm for UAV-mounted base stations,” in _IEEE Wireless Communications and Networking Conference (WCNC 2019)_ , 2019, pp. 1–6.
  * [14] N. Nouri, J. Abouei, A. R. Sepasian, M. Jaseemuddin, A. Anpalagan, and K. N. Plataniotis, “Three-dimensional multi-UAV placement and resource allocation for energy-efficient IoT communication,” _IEEE Internet of Things Journal_ , vol. 9, no. 3, pp. 2134–2152, 2022.
  * [15] R. Chen, Y. Sun, L. Liang, and W. Cheng, “Joint power allocation and placement scheme for UAV-assisted IoT with QoS guarantee,” _IEEE Transactions on Vehicular Technology_ , vol. 71, no. 1, pp. 1066–1071, 2022.
  * [16] M. Alzenad, A. El-Keyi, F. Lagum, and H. Yanikomeroglu, “3-D placement of an unmanned aerial vehicle base station (UAV-BS) for energy-efficient maximal coverage,” _IEEE Wireless Communications Letters_ , vol. 6, no. 4, pp. 434–437, 2017.
  * [17] V. Saxena, J. Jaldén, and H. Klessig, “Optimal UAV base station trajectories using flow-level models for reinforcement learning,” _IEEE Transactions on Cognitive Communications and Networking_ , vol. 5, no. 4, pp. 1101–1112, 2019.
  * [18] A. Al-Hourani, S. Kandeepan, and S. Lardner, “Optimal LAP altitude for maximum coverage,” _IEEE Wireless Communications Letters_ , vol. 3, no. 6, pp. 569–572, 2014.
  * [19] L.-Y. Wu, X.-S. Zhang, and J.-L. Zhang, “Capacitated facility location problem with general setup cost,” _Computers & Operations Research_, vol. 33, no. 5, 2006.
  * [20] R. Saber and R. Murray, “Flocking with obstacle avoidance: cooperation with limited communication in mobile networks,” in _42nd IEEE International Conference on Decision and Control (CDC 2003)_ , Maui, HI, USA 2003.
  * [21] R. Olfati-Saber, “Flocking for multi-agent dynamic systems: Algorithms and theory,” _IEEE Transactions on Automatic Control_ , vol. 51, no. 3, 2006\.
  * [22] M. J. Farooq and Q. Zhu, “A multi-layer feedback system approach to resilient connectivity of remotely deployed mobile Internet of things,” _IEEE Transactions on Cognitive Communications and Networking_ , vol. 4, no. 2, pp. 422–432, 2018.
  * [23] ——, “Cognitive connectivity resilience in multi-layer remotely deployed mobile internet of things,” in _GLOBECOM 2017 - 2017 IEEE Global Communications Conference_ , 2017, pp. 1–6.


11institutetext: Purple Mountain Observatory and Key Laboratory of Radio
Astronomy, Chinese Academy of Sciences, 10 Yuanhua Road, Nanjing 210023, China
11email<EMAIL_ADDRESS>22institutetext: School of Astronomy & Space
Science, Nanjing University, 163 Xianlin Avenue, Nanjing 210023, China
33institutetext: Key Laboratory of Modern Astronomy and Astrophysics, Nanjing
University, Ministry of Education, Nanjing 210023, China

# Discovery of an old supernova remnant candidate through carbon monoxide line
emission

Xin Zhou 11 Yang Su 11 Ji Yang 11 Yang Chen 2233 Zhibo Jiang 11

Most old supernova remnants (SNRs) in the Milky Way have not yet been
identified. Considering their large potential number and the sufficient
momentum-energy transfer to the interstellar medium (ISM), they are a key part
of our understanding of the overall role of SNRs in the ISM. Here we report
our discovery of an expanding molecular shell identified by CO line
observations, namely G16.11$-$0.51. It covers a known SNR, specifically
G16.0$-$0.5, and is larger in size, i.e. 0.56∘ over 0.20∘. Based on its
spatial and kinematic structures, weak nonthermal radio-continuum emission,
and derived physical properties, we suggest that it is an old SNR. At a
systemic velocity of +41.3 km s-1, the best estimated kinematic distance of
G16.11$-$0.51 is $\sim$3.2 kpc, implying its radius of about 15.6 pc. The age
of G16.11$-$0.51 is estimated to be greater than $\sim$$10^{5}$ yr, and, in a
dense molecular environment, it has formed dense and thin shell layers. The
kinetic energy of the expanding molecular gas of G16.11$-$0.51 is about
$6.4\times 10^{49}$ erg, accounting for approximately six percent of the
initial SN explosion energy. Although old SNRs have essentially become cold
and hard to detect, our discovery suggests that they can be found by searching
for CO line emissions.

###### Key Words.:

ISM: supernova remnants – ISM: molecules– ISM: clouds – ISM: lines and bands

##  1 Introduction

Supernova remnants (SNRs) contain a large amount of momentum, energy, and
heavy elements. They sufficiently transfer them to the surrounding
interstellar medium (ISM) and cool down in their late stages (e.g., $>\sim
10^{5}$; Koo et al., 2020). One supernova event occurs on average every $40\pm
10$ yr in the Milky Way (Tammann et al., 1994), resulting in over 2000 SNRs
younger than $10^{5}$ yr. If we take $10^{6}$ yr to be the lifetime of an SNR,
there would be over 20000 SNRs in total. The number of old SNRs is much larger
than that of young bright ones. Old SNRs are a key part of understanding the
overall role of SNRs in the ISM, especially their impact on molecular clouds
(MCs; e.g., how they regulate the formation of next generation stars). Old
SNRs that have a large potential number and large size are also likely to
overlap with other sources, and they act as backgrounds to affect our
understanding of other sources.

Only a limited number of SNRs have been identified in our Galaxy, i.e. fewer
than 400, mostly based on their radio-continuum emission (cf. Green, 2019;
Ferrand & Safi-Harb, 2012). Most of these known SNRs are bright and young. It
is difficult to detect old SNRs because they usually have very weak radiation;
for example, they are radio faint and begin to dissolve into the ISM.
Benefiting from enhanced observations in different wavelengths, new SNR
candidates have been continuously discovered, e.g., by optical emission (Fesen
et al., 2020), by X-ray emission (Becker et al., 2021; Churazov et al., 2021;
Khabibullin et al., 2023), and mostly by radio-continuum emission (Anderson et
al., 2017; Hurley-Walker et al., 2019a; Gao et al., 2020; Dokara et al., 2021;
Ball et al., 2023, etc.). It appears that there are still a large number of
potential SNRs waiting to be discovered; nevertheless, background emission
would limit such detections (see, for example, radio-continuum emission from
HESS J1912+101 confused by background diffuse emission Reich & Sun, 2019),
especially for weak SNRs or toward the inner Galactic region. Additionally,
the large size and quantity of old SNRs can lead to a severe overlapping
effect, further complicating the situation. Considering that old SNRs would
form a cold, dense, and gradually momentum-conserving shell, H i 21 cm line
emission has been used as a good tracer to detect such old SNRs, especially
large ones at high Galactic latitude or with localized high-velocity H i
features (e.g., Koo & Kang, 2004; Koo et al., 2006; Kang & Koo, 2007; Kang et
al., 2012; Xiao & Zhu, 2014, etc.). Although CO line emission can also trace
dense shells associated with SNRs (e.g., for known SNRs; Chen et al., 2014;
Sofue et al., 2021), no known SNR has been discovered through CO lines to date
(Green, 2019; Ferrand & Safi-Harb, 2012). Note that, based on the analysis of
CO and H i data, Su et al. (2017) suggested the SNR origin of the TeV source
HESS J1912+101. Compared to H i gas, molecular gases are more dense and
clumpy, hence, it is difficult to detect a large complete shell-like structure
associated with an SNR. In addition, molecular lines are very likely to have
serious overlapping effects too, especially toward the inner Galactic region,
which can complicate their velocity features. Nevertheless, since remnants of
core-collapse supernovae are thought to be located close to their parent MCs,
numerous molecular shells are expected to have originated from old SNRs, e.g.,
some CO bubbles in the W43 molecular complex were interpreted as fully evolved
SNRs by Sofue (2021).

Here we report the discovery of a new expanding CO shell, namely
G16.11$-$0.51, which covers a known SNR, i.e. G16.0$-$0.5. Its spatial and
kinematic structures, weak nonthermal radio-continuum emission, and derived
physical properties indicate its origin from an SN explosion.

##  2 Observations

The CO line emission was observed during November 2011 to November 2013 and
from February to April 2019 using the Purple Mountain Observatory (PMO) 13.7 m
millimeter-wavelength telescope located in Delingha, China, as part of the
Milky Way Imaging Scroll Painting (MWISP)
project111http://english.dlh.pmo.cas.cn/ic/. The $3\times 3$ multibeam
sideband separation receiver, i.e. Superconducting Spectroscopic Array
Receiver (SSAR, Shan et al., 2012), was used to simultaneously observe the
12CO (J=1–0), 13CO (J=1–0), and C18O (J=1–0) lines. The fast Fourier transform
(FFT) spectrometer with 1 GHz bandwidth and 16384 channels was used as the
backend for each sideband. Correspondingly, the spectral resolutions of the
three CO lines were 0.17 km s-1 for 12CO (J=1–0) and 0.16 km s-1 for both 13CO
(J=1–0) and C18O (J=1–0). We mapped a $1.2^{\circ}\times 1.2^{\circ}$ area via
on-the-fly (OTF) observing mode. Following the standard procedure of the MWISP
CO line survey, observations were performed at least twice along the Galactic
longitude and latitude directions, respectively, and specialized data checks
were carried out. Data with bad baselines were excluded, and thereby baseline
problems were largely suppressed. However, some minor rms noise
inhomogeneities were introduced, for example, along the Galactic longitude or
latitude direction. The total error in pointing and tracking was within
$5^{\prime\prime}$. The data were meshed with a grid spacing of
$30^{\prime\prime}$., and the half-power beam width (HPBW) was $\sim 51"$. The
typical main-beam efficiencies $\eta_{\rm mb}$ were $\sim$48% for USB and
$\sim$52% for LSB, which was used to convert the antenna temperature to the
main beam temperature by $T_{\rm mb}=T_{A}^{*}/\eta_{\rm mb}$ (see details in
the status report of the telescope222http://www.radioast.nsdc.cn/mwisp.php).
After performing the linear baseline subtraction process, the typical rms
noises obtained were $\sim 0.43$ K per channel for 12CO (J=1–0) and $\sim
0.21$ K per channel for 13CO (J=1–0) and C18O (J=1–0). All data were reduced
using the GILDAS/CLASS package333http://www.iram.fr/IRAMFR/GILDAS.

Radio continuum emission data in four bands with frequencies ranging from 72
to 103 MHz, 103 to 134 MHz, 138 to 170 MHz, and 170 to 231 MHz were obtained
from the Galactic and Extragalactic All-sky Murchison Widefield Array (GLEAM)
survey444https://www.mwatelescope.org/gleam (Wayth et al., 2015; Hurley-Walker
et al., 2017, 2019b). Other radio-continuum emission data in six bands were
obtained from The HI/OH/Recombination line survey (THOR; Beuther et al., 2016;
Wang et al., 2020), which are centered at 1.06, 1.31, 1.44, 1.69, 1.82, and
1.95 GHz and have a bandwidth of 128 MHz. Combined 1.4 GHz continuum data from
the THOR and the VLA Galactic Plane Survey (VGPS, Stil et al., 2006) was also
applied, from which the flux retrieved is consistent with the literature
(Anderson et al., 2017; Wang et al., 2020).

##  3 Results

Figure 1: Integrated intensity and position-velocity maps of CO emission
toward G16.11$-$0.51. Bottom left panel: pseudo tricolor image of the
integrated intensity of 12CO (J=1–0) (blue), 13CO (J=1–0) (green), and
C18O(J=1–0) (red) line emission in the velocity range from +35 to +85 km s-1.
The intensities of all three CO lines have confidence levels over $5\sigma$.
The large and small yellow dashed circles indicate the extent of G16.11$-$0.51
and SNR G16.0$-$0.5, respectively. The beam is represented by a white circle
in the lower left corner. The top and right panels show the 12CO (J=1–0)
position-velocity maps along the strips indicated in the intensity map. The
minimum value of the background grayscale map is $1\sigma$, and the contour
levels are from $3\sigma$ and in a step of $1\sigma$. Black ellipses indicate
structures of expanding molecular gas around the boundary of G16.11$-$0.51.
Three red arrows in the right panel mark the locations of three examples of
broadened CO lines associated with G16.11$-$0.51.

We show the pseudo tricolor image of integrated intensities of three CO
isotope lines of G16.11$-$0.51 over a wide velocity range in Figure 1.
Multiple shell structures are clearly seen around a circular region with a
radius of about $0.28^{\circ}$, which are distributed from the north to the
southwest and in the east. G16.11$-$0.51 overlaps with a known SNR G16.0$-$0.5
but it is much larger in size. The molecular shell associated with G16.0$-$0.5
in its southern region is also visible, which was detected in previous work
too (Beaumont et al., 2011). As shown in the position-velocity maps along
horizontal and vertical strips across G16.11$-$0.51 (Figure 1), the molecular
shells are expanding at a velocity of about 10 km s-1. Some molecular clumps
correlated with the expanding shells are also strongly disturbed, indicated by
line broadening features of about 20 km s-1 (see, for example, those indicated
by three red arrows in the right panel of Figure 1). We note that the line
broadening features are significant and not due to scanning effects (artifacts
caused by baseline variations). Two of these molecular clumps are located
within and on the edge of the G16.0$-$0.5 region. However, the blueshifted
velocity of these molecular clumps is consistent with the blueshifted velocity
direction of the expanding shell at their locations. The quiescent molecular
gas at $\sim$+43 km s-1 around the expanding shell structures is rarefied.
Molecular gases at $\sim$+57 km s-1 inside G16.11$-$0.51 are all disturbed or
associated with the expanding shells. There is a small amount of quiescent
molecular gas at $\sim$+57 km s-1 outside G16.11$-$0.51. However, this small
amount of molecular gas may have a peculiar velocity compared to the Galactic
rotation curve, but the major $\sim$+43 km s-1 component is unlikely. The
expanding shell structures and broadened line emission are probably associated
with the quiescent molecular gas at $\sim$+43 km s-1. The $\sim$+43 and
$\sim$+57 km s-1 MCs may be adjacent to each other and at the same distance,
however, it is uncertain.

Figure 2:  Velocity dispersion (moment 2) maps of 12CO (left) and 13CO (right)
line emission over the velocity range from +40 to +60 km s-1, overlaid with
the GLEAM 170–230 MHz radio-continuum emission contours
($\sim$2${}^{\prime}.25$ angular resolution). The contour levels are 0.76 and
0.81 Jy/beam. The black dashed circles are the same as in Figure 1, indicating
the extent of G16.11$-$0.51 and SNR G16.0$-$0.5, respectively. The H ii
regions are indicated by magenta circles in the right panel, and all are
introduced from the WISE catalog of Galactic H ii regions (Anderson et al.,
2014), which is one of the most complete catalogs of H ii regions in the
Galaxy. Radio continuum sources extracted from the THOR survey are marked by
black solid boxes (Bihr et al., 2016; Wang et al., 2018).

The kinematic properties of molecular gases and the corresponding spatial
distribution can be illustrated by the velocity dispersion map. As shown in
the 12CO velocity dispersion map (left panel, Figure 2), disturbed molecular
gases are present in the northwestern shell of G16.11$-$0.51 and in a
molecular clump located around the southeastern boundary. There is also some
diffuse and disturbed molecular gas distributed around the northeastern
boundary of G16.11$-$0.51, where there is a weak radio-continuum shell. Known
H ii regions are introduced from the WISE catalog of Galactic H ii regions
(Anderson et al., 2014), which is one of the most complete catalogs of H ii
regions in the Galaxy (see the right panel of Figure 2). Note that the
systemic velocity of some H ii regions is not within the selected velocity
range. The disturbed molecular gases along the G16.11$-$0.51 boundary do not
originate from known H ii regions. For SNR G16.0$-$0.5, disturbed gas in the
associated molecular shell is well correlated with the southeastern radio-
continuum shell of the remnant. The 13CO line emission is optically thinner
than the 12CO line emission, which can trace inner and denser molecular gas.
The 13CO velocity dispersion map (right panel, Figure 2) shows the
distribution of dense disturbed molecular gas. By comparison, no dense
disturbed gas associated with G16.0$-$0.5 is seen. For G16.11$-$0.51, dense
disturbed gas distributed in a thin layer is seen inside its northwestern
molecular shell. Besides this, in its southeastern clump, dense disturbed gas
is distributed only in the thin inner layer facing the center of
G16.11$-$0.51. It indicates that the disturbed gas in the southeastern clump
is mainly distributed on its surface and transmitted toward the outside of
G16.11$-$0.51. Further study of these disturbed molecular gases in future work
is needed for a better understanding of them. For MCs outside the
G16.11$-$0.51 region, some disturbed gases are associated with H ii regions,
such as those in the northwest, northeast, and southeast.

Figure 3:  Flux densities of G16.11$-$0.51 obtained from the GLEAM data only
(left) and the GLEAM+THOR+(THOR+VGPS) data (right). The flux density obtained
from the THOR+VGPS 1.4 GHz continuum data is presented by a blue star in the
right panel. The flux densities are extracted from the G16.11$-$0.51 region
indicated in Figure 2. An annular region half the beamsize of the GLEAM 72–103
MHz data away from the source region, with an area similar to the source
region, is used to subtract the background emission. The flux densities of H
ii regions, THOR radio-continuum sources, and the SNR G16.0$-$0.5 region, all
enlarged by half the beamsize, were removed from the calculation (see regions
in the right panel of Figure 2). The axes in the left panel are logarithmic,
while the ones in the right panel are linear. Internal and external flux
density scale errors are applied as $2\%$ and $8\%$ for the GLEAM data in left
and right panels, respectively (Hurley-Walker et al., 2017). The flux density
uncertainty of the THOR and THOR+VGPS data is determined by measuring the
variation of the emission free region (see Anderson et al., 2017, for
reference). The fit power-law spectra are shown by solid lines. The best-fit
spectral indices are $-0.8\pm 0.3$ and $-0.77\pm 0.04$ for the GLEAM and
GLEAM+(THOR+VGPS) data, respectively. The upper limits of the flux densities
from the THOR-only data are shown for comparison but are not used in fitting.

G16.11$-$0.51 presents a partial radio-continuum shell along its northeastern
boundary (see Figure 2). There is significant radio-continuum emission outside
G16.11$-$0.51 in the northeast, which is associated with H ii regions (see the
right panel of Figure 2). However, radio-continuum emissions from
G16.11$-$0.51 and the H ii regions surround different centers. We measured the
flux densities of G16.11$-$0.51 in different bands, with the background level
subtracted. We first smoothed all the radio-continuum data in different bands
to a common resolution as the lowest resolution GLEAM 72–103 MHz data. An
annular region half the beamsize of the GLEAM 72–103 MHz data away from the
source region, with an area similar to the source region, was used to subtract
the background emission. The flux densities of H ii regions, THOR radio-
continuum sources, and the SNR G16.0$-$0.5 region, all enlarged by half the
beamsize, are removed from the calculation (see regions in the right panel of
Figure 2). The same regions are applied for all GLEAM, THOR, and THOR+VGPS
bands. The radio-continuum emission of G16.11$-$0.51 is weak, hence, the
individual measured flux densities are not so significant. The overall radio-
continuum spectrum can be fit by the power-law function (see Figure 3). The
best-fit spectral indices obtained using the GLEAM-only and GLEAM+(THOR+VGPS)
data are consistent and steep, i.e. $\alpha=-0.8\pm 0.3$
($S_{\nu}\propto\nu^{\alpha}$). No radio-continuum emission is detected by the
THOR-only data, which is likely to suffer from serious missing flux problems
in the interferometric observations for such a large source. These indicate
that the radio-continuum emission is nonthermal and not from H ii regions. The
extrapolated 1 GHz surface brightness of G16.11$-$0.51 is estimated as
$\Sigma_{1\leavevmode\nobreak\ GHz}\sim 3\times 10^{-21}$ W m-2 Hz-1 sr-1.

##  4 Discussion

Figure 4: 12CO (J=1–0) (blue) and 13CO (J=1–0) (green) spectra extracted from
1.2 times enlarged G16.11$-$0.51 region with 1.2 times enlarged SNR
G16.0$-$0.5 region subtracted, together with their best-fit Gaussian model and
residuals. The 12CO components peak at +19.3, +29.5, +43.3, +58.0, +68.1, and
+76.0 km s-1, and the 13CO components peak at +19.5, +30.2, +41.3, and +57.4
km s-1. Individual Gaussian models of the $\sim$+43 and $\sim$+57 km s-1
components are shown by black and red dotted lines, respectively. $5\sigma$
levels of the residuals are shown by blue and green dotted lines for 12CO and
13CO, respectively. The 13CO (J=1–0) spectrum and its fitting result and
residuals are multiplied by a factor of 2 for better visibility.  Table 1: Fit
and derived parameters for two related velocity components in G16.11$-$0.51.

Gaussian Components |
---|---
Component | Line |  Peak $T_{\rm mb}$ |  Center $V_{\rm LSR}$aa$a$$V_{\rm LSR}$ is the velocity with respect to the local standard of rest. | FWHM |
|  | (K) |  ( km s-1) |  (km s-1) |
rest |  12CO (J=1–0) | $4.79\pm 0.03$ | $+43.31\pm 0.07$ | $12.4\pm 0.2$ |
|  13CO (J=1–0) | $1.027\pm 0.007$ | $+41.25\pm 0.09$ | $12.0\pm 0.3$ |
shifted |  12CO (J=1–0) | $1.30\pm 0.03$ | $+58.0\pm 0.2$ | $8.8\pm 0.4$ |
|  13CO (J=1–0) | $0.198\pm 0.008$ | $+57.4\pm 0.2$ | $7.8\pm 0.5$ |
Derived Physical Parameters
Component |  $T_{\rm ex}$bb$b$Using the assumption of local thermal equilibrium (LTE). See the details of the calculation method in Zhou et al. (2016). |  $\tau(^{13}$CO)bb$b$Using the assumption of local thermal equilibrium (LTE). See the details of the calculation method in Zhou et al. (2016). |  $N({\rm H_{2}})$cc$c$Derived from 13CO column density by assuming the 13CO abundance of 1.4$\times 10^{-6}$ (Ripple et al., 2013). For comparison, we also show the values in brackets, which are estimated by using the conversion factor $N({\rm H}_{2})/{W}(^{12}{\rm CO)}\simeq 1.8\times 10^{20}\leavevmode\nobreak\ {\rm cm}^{-2}\leavevmode\nobreak\ {\rm K}^{-1}\leavevmode\nobreak\ {\rm km}^{-1}$ s (Dame et al., 2001). $d_{3.2}$ stands for $d/(3.2\leavevmode\nobreak\ {\rm kpc})$, where d is the distance to G16.11$-$0.51 in units of kpc. |  $M$cc$c$Derived from 13CO column density by assuming the 13CO abundance of 1.4$\times 10^{-6}$ (Ripple et al., 2013). For comparison, we also show the values in brackets, which are estimated by using the conversion factor $N({\rm H}_{2})/{W}(^{12}{\rm CO)}\simeq 1.8\times 10^{20}\leavevmode\nobreak\ {\rm cm}^{-2}\leavevmode\nobreak\ {\rm K}^{-1}\leavevmode\nobreak\ {\rm km}^{-1}$ s (Dame et al., 2001). $d_{3.2}$ stands for $d/(3.2\leavevmode\nobreak\ {\rm kpc})$, where d is the distance to G16.11$-$0.51 in units of kpc. |
| (K) |  | $(10^{21}\leavevmode\nobreak\ \,{\rm cm}^{-2})$ | $(10^{4}\leavevmode\nobreak\ {\rm M}_{\rm\odot})$ |
rest | $8.1$ | $0.27$ |  $11$ (11) |  $14$$d_{3.2}$2 ($17$$d_{3.2}$2) |
shifted | $4.3$ | $0.19$ |  $1.9$ (2.2) |  $2.5$$d_{3.2}$2 ($3.2$$d_{3.2}$2) |
555

We extracted CO spectra in the G16.11$-$0.51 region, which is defined as a
circular region 1.2 times larger than that shown in Figure 2. We also excluded
the SNR G16.0$-$0.5 region, which is a circular region 1.2 times larger than
that shown in Figure 2. The extracted 12CO and 13CO spectra are shown in
Figure 4. According to the number of their peaks, the 12CO spectrum is fit by
six Gaussian functions, and the 13CO spectrum is fit by four Gaussian
functions. As shown in Section 3, G16.11$-$0.51 is associated with the MC at a
systemic velocity of +41.3 km s-1. Molecular gas at +57.4 km s-1 within the
G16.11$-$0.51 region is mainly driven by G16.11$-$0.51 from the +41.3 km s-1
MC. The different velocities between the expanding shell of G16.11$-$0.51 and
its parent MC indicate that the medium surrounding G16.11$-$0.51 has a density
gradient in the line of sight. The schematic diagram of the formation of the
expanding shell of G16.11$-$0.51 is shown in the appendix. The fit and derived
parameters for the +41.3 and +57.4 km s-1 components are listed in Table 1.
The kinematic distances of the +41.3 km s-1 MC are estimated as $3.2\pm 0.4$
and $4.0\pm 0.3$ kpc, based on a full distance probability density
function666http://bessel.vlbi-astrometry.org/node/378 (Reid et al., 2016,
2019). The +41.3 km s-1 MC is more likely to be at 3.2 kpc, which has a
probability of 0.53. The momentum and kinetic energy of the expanding
molecular gas of G16.11$-$0.51 are estimated to be $4\times
10^{5}$$d_{3.2}$${}^{2}\leavevmode\nobreak\ {\rm M}_{\rm\odot}\,{\rm km}\,{\rm
s}^{-1}$ and $6.4\times 10^{49}$$d_{3.2}$2 erg, respectively. We note that,
although a far kinematic distance of the +41.3 km s-1 MC can be estimated to
be $\sim$12 kpc, the corresponding kinetic energy of the molecular gas would
be too large, i.e. $9\times 10^{50}$ erg, even comparable to the typical SN
explosion energy value. If the expanding shell is produced by the stellar wind
of the OB star, even if all the stellar wind energy is transferred to the
shell, a star with a mass greater than that of an O9 V type star is required
(e.g., Abbott, 1982). However, we searched the SIMBAD astronomical database
(Wenger et al., 2000) within a 1.2 times enlarged region of G16.11$-$0.51 and
found no O-type stars. Only a few B-type stars are found in the region, which
are located at distances of less than $\sim$1.2 kpc, and, hence, they are not
related to G16.11$-$0.51. The expanding shell with thin and dense layers, and
its physical properties, indicates an SN origin for G16.11$-$0.51, which is
further supported by the associated nonthermal radio-continuum emission.

If G16.11$-$0.51 is at a distance of 3.2 kpc, its radius is 15.6 pc. The weak
radio-continuum emission as well as the presence of dense shell structures
indicate that SNR G16.11$-$0.51 is in the radiative stage. Since the molecular
shells are basically distributed on a circle on the plane of the sky, some
simple models of the evolution of SNRs in a homogeneous medium can be applied
as approximations to provide some references. Assuming that the SNR evolved in
a homogeneous ISM, its age can be estimated as $t=2r_{s}/(7v_{s})\sim
2.7\times 10^{5}$$d_{3.2}$ yr (McKee & Ostriker, 1977), and the explosion
energy as $E_{SN}=6.8\times 10^{43}n_{0}^{1.16}(v_{s}/1\leavevmode\nobreak\
{\rm km\leavevmode\nobreak\ s}^{-1})^{1.35}(r_{s}/1\leavevmode\nobreak\ {\rm
pc})^{3.16}\zeta_{m}^{0.161}\sim 2\times 10^{51}$$d_{3.2}$2 erg, where the
ambient hydrogen density $n_{0}$ is estimated to be $\sim$63 cm-3, the shock
velocity $v_{s}$ as $v_{mol}=16.1$ km s-1, the radius of the SNR $r_{s}$ 15.6
pc, and the metallicity parameter $\zeta_{m}=Z/Z_{\rm\odot}$ is set to be 1
(Cioffi et al., 1988). The ambient density is estimated by assuming that the
mass of the expanding molecular shell was initially uniformly distributed
throughout its internal volume. The explosion energy obtained is big in
comparison to the typical value, i.e. $10^{51}$ erg. As the SNR evolves in a
non-uniform medium, there is a leakage of the SNR energy in less dense
directions, and therefore the explosion energy is even underestimated.
Considering that the SNR is likely to have evolved in a bubble blown by the
stellar wind of its progenitor, with the contribution of the stellar wind to
the accumulation of the ISM, less SN explosion energy will be required. In the
wind-blown bubble scenario, we estimate the explosion energy and age using the
model of Chen et al. (2003). Assuming $\lambda=r_{j}/r_{s}=0.9$ with
$(\eta,\beta)=(1,0)$, its explosion energy can be estimated as
$E_{SN}=1.05\rho_{0}v_{s}^{2}r_{s}^{3}\lambda^{-2}[F_{v}^{(R)}(\lambda)]^{-2}\sim
10^{51}$$d_{3.2}$2 erg and the age as
$t=1.02\rho_{0}^{0.5}r_{s}^{2.5}E_{SN}^{-0.5}\lambda^{-1}F_{r}^{(R)}(\lambda)\sim
10^{5}$$d_{3.2}$ yr, where $r_{j}$ is the radius of the progenitor’s wind-
blown bubble, $\rho_{0}$ is the density of the wall of the bubble, the
dimensionless term $F_{v}^{(R)}(\lambda)$ is estimated to be $\sim$0.444, and
$F_{r}^{(R)}(\lambda)$ is estimated to be $\sim$0.206 (cf. Equations [26],
[29], [31], and [35] in Chen et al., 2003). For $\lambda=r_{j}/r_{s}=0.9$, the
radius of the wind-blown bubble is $r_{j}\sim 14$ pc. The assumption of
$\eta=1$ means that the SNR enters the radiative phase directly after hitting
the cavity wall, which is reasonable for $r_{j}\gtrsim 3.8$ pc (see Section
2.1 in Chen et al., 2003), and the assumption of $\beta=0$ means that the
density inside the cavity is much smaller than that of the cavity wall. The
density of the wall of the wind-blown bubble is estimated by assuming that the
mass of the expanding molecular shell was initially evenly distributed in a
spherical shell with inner and outer radii of $r_{j}$ and $r_{s}$,
respectively, i.e. $n_{0}$$\sim$233 cm-3. According to the radius of the wind-
blown bubble, the progenitor’s initial mass can be estimated to be $\sim$19
${\rm M}_{\rm\odot}$ by applying the linear relationship between the size of a
massive star’s main-sequence bubble in a molecular environment and its initial
mass (Chen et al., 2013). In the wind-blown bubble scenario, a relatively
reasonable explosion energy can be obtained. Note that the uniform wind-blown
bubble model is still a rough approximation. Since the radius of the SNR on
the plane of the sky is smaller than the average radius (see the appendix for
more information), the age obtained (i.e., $\sim$$10^{5}$ yr) is only a lower
limit. Further analysis in future work is needed to obtain more appropriate
models based on specific molecular gas distributions (e.g., numerical models).

The radio-continuum emission from G16.11$-$0.51 is weak, of which the diameter
($\sim$31 pc) and the radio-continuum surface brightness
($\Sigma_{1\leavevmode\nobreak\ GHz}\sim 3\times 10^{-21}$ W m-2 Hz-1 sr-1)
conform to the $\Sigma$–D relation (Pavlović et al., 2018). In general, as the
SNRs become older, their radio-continuum emissions become weaker, e.g., the
SNRs discovered by H i observations with an age of around one million years is
with $\Sigma_{1GHz}\lesssim 7\times 10^{-23}$ W m-2 Hz-1 sr-1(Koo et al.,
2006; Xiao & Zhu, 2014). Nevertheless, when the SNR shock encounters dense
clouds, due to the compressed magnetic fields and the accelerated preexisting
cosmic ray electron, the radio-continuum emission can be enhanced (Blandford &
Cowie, 1982; Draine & McKee, 1993). It seems to be efficient for SNRs
encountering clouds while still having quite a high Mach number (Pavlović et
al., 2018). No radio-continuum emission is detected around the northwestern
and southeastern molecular shells, which are probably swept up by the SNR
shock but not preexisting ones. The nonthermal spectral index of the radio-
continuum emission from G16.11$-$0.51 is steep, which resembles indices of
shell-type SNRs detected in the GLEAM survey (e.g., Hurley-Walker et al.,
2019a). The steep spectral index can be expected for an old SNR due to the
energy loss of electrons through synchrotron radiation (e.g., Xiao & Zhu,
2014).

G16.11$-$0.51 completely covers SNR G16.0$-$0.5 and is larger in size.
Beaumont et al. (2011) used a machine learning algorithm to classify broadened
12CO (J=3–2) line structures in a position-position-velocity data cube for SNR
G16.0$-$0.5. Most of the coverage of their data is in the G16.11$-$0.51
region. In addition to most of the identified broad lines distributed within
SNR G16.0$-$0.5, there are also broad lines identified outside G16.0$-$0.5 but
within G16.11$-$0.51 (see, for example, Figure 8 in Beaumont et al., 2011). As
shown in Section 3 (Figure 1) here, there are broad lines both inside and
outside the G16.0$-$0.5 region correlated with the expanding CO shell of
G16.11$-$0.51. Even if it is correlated with G16.11$-$0.51, we cannot totally
rule out the possibility that the broad lines in the G16.0$-$0.5 region are
associated with SNR G16.0$-$0.5. For old SNRs, their shock waves are weak, and
it is difficult to accelerate the molecular gas to very high velocities. The
molecular clumps with broad line features, which are correlated with
G16.11$-$0.51, may have been accelerated earlier by the fast SNR shock and
they may have experienced less deceleration due to the greater density. The
molecular gas associated with SNR G16.0$-$0.5 is mainly at $\sim$+68 km s-1,
and the $\sim$+68 km s-1 velocity component is an independent MC distributed
both inside and outside SNR G16.0$-$0.5 (see Zhou et al., 2023). Accordingly,
most of the broadened CO lines associated with SNR G16.0$-$0.5 are
blueshifted. If the molecular gas at $\sim$+68 km s-1 within SNR G16.0$-$0.5
is redshifted from the +41.3 km s-1 component by the stellar wind of its
progenitor, it should be further redshifted by the shock wave of the
subsequent remnant, which is not the case. It is still possible that SNR
G16.0$-$0.5 is located near the shell on the farside of G16.11$-$0.51, but
there is no evidence to support this.

In conclusion, we suggest G16.11$-$0.51 is an old SNR. As a background to SNR
G16.0$-$0.5, some broad lines in the G16.0$-$0.5 region show correlations with
G16.11$-$0.51. The finding of such an SNR characterized by an expanding CO
shell may suggest other similar remnants waiting to be discovered.

###### Acknowledgements.

We thank Ping Zhou for helpful feedback on the draft manuscript which improved
the paper. We also thank the anonymous referee for very valuable comments that
improved this paper and its conclusions. We are grateful to all the members of
the Milky Way Scroll Painting-CO line survey group, especially the staff of
the Qinghai Radio Observing Station at Delingha, for their support during the
observation. In particular, we thank the anonymous referee of Zhou et al.
(2023) for providing very helpful suggestions that inspired our in-depth study
of this region. This work is part of the Milky Way Imaging Scroll Painting
(MWISP) multi-line survey project, which is supported by the National Key R&D
Program of China grant no. 2023YFA1608000 and 2017YFA0402701, and the Key
Research Program of Frontier Sciences, CAS, grant no. QYZDJ-SSW-SLH047. Y.S.,
J.Y., and Y.C. acknowledge support from the NSFC grants 12173090, 12041305,
12173018 and 12121003. This research has made use of the SIMBAD database
operated at CDS, Strasbourg, France. We acknowledge the use of the VGPS data;
the National Radio Astronomy Observatory is a facility of the National Science
Foundation operated under cooperative agreement by Associated Universities,
Inc.

### Data Availability

The MWISP CO (J=1–0) data of SNR G16.11$-$0.51 is available online in
https://www.scidb.cn/en, at https://doi.org/10.57760/sciencedb.15577.

## References

  * Abbott (1982) Abbott, D. C. 1982, ApJ, 263, 723
  * Ambrocio-Cruz et al. (2006) Ambrocio-Cruz, P., Rosado, M., & de La Fuente, E. 2006, Rev. Mexicana Astron. Astrofis., 42, 241
  * Anderson et al. (2014) Anderson, L. D., Bania, T. M., Balser, D. S., et al. 2014, ApJS, 212, 1
  * Anderson et al. (2017) Anderson, L. D., Wang, Y., Bihr, S., et al. 2017, A&A, 605, A58
  * Ball et al. (2023) Ball, B. D., Kothes, R., Rosolowsky, E., et al. 2023, MNRAS, 524, 1396
  * Beaumont et al. (2011) Beaumont, C. N., Williams, J. P., & Goodman, A. A. 2011, ApJ, 741, 14
  * Becker et al. (2021) Becker, W., Hurley-Walker, N., Weinberger, C., et al. 2021, A&A, 648, A30
  * Beuther et al. (2016) Beuther, H., Bihr, S., Rugel, M., et al. 2016, A&A, 595, A32
  * Bihr et al. (2016) Bihr, S., Johnston, K. G., Beuther, H., et al. 2016, A&A, 588, A97
  * Blandford & Cowie (1982) Blandford, R. D. & Cowie, L. L. 1982, ApJ, 260, 625
  * Chen et al. (2014) Chen, Y., Jiang, B., Zhou, P., et al. 2014, in IAU Symposium, Vol. 296, Supernova Environmental Impacts, ed. A. Ray & R. A. McCray, 170–177
  * Chen et al. (2003) Chen, Y., Zhang, F., Williams, R. M., & Wang, Q. D. 2003, ApJ, 595, 227
  * Chen et al. (2013) Chen, Y., Zhou, P., & Chu, Y.-H. 2013, ApJ, 769, L16
  * Churazov et al. (2021) Churazov, E. M., Khabibullin, I. I., Bykov, A. M., et al. 2021, MNRAS, 507, 971
  * Cioffi et al. (1988) Cioffi, D. F., McKee, C. F., & Bertschinger, E. 1988, ApJ, 334, 252
  * Dame et al. (2001) Dame, T. M., Hartmann, D., & Thaddeus, P. 2001, ApJ, 547, 792
  * Dokara et al. (2021) Dokara, R., Brunthaler, A., Menten, K. M., et al. 2021, A&A, 651, A86
  * Draine & McKee (1993) Draine, B. T. & McKee, C. F. 1993, ARA&A, 31, 373
  * Ferrand & Safi-Harb (2012) Ferrand, G. & Safi-Harb, S. 2012, Advances in Space Research, 49, 1313
  * Fesen et al. (2020) Fesen, R. A., Weil, K. E., Raymond, J. C., et al. 2020, MNRAS, 498, 5194
  * Gao et al. (2020) Gao, X. Y., Reich, P., Reich, W., Hou, L. G., & Han, J. L. 2020, MNRAS, 493, 2188
  * Green (2019) Green, D. A. 2019, Journal of Astrophysics and Astronomy, 40, 36
  * Hurley-Walker et al. (2017) Hurley-Walker, N., Callingham, J. R., Hancock, P. J., et al. 2017, MNRAS, 464, 1146
  * Hurley-Walker et al. (2019a) Hurley-Walker, N., Filipović, M. D., Gaensler, B. M., et al. 2019a, PASA, 36, e045
  * Hurley-Walker et al. (2019b) Hurley-Walker, N., Hancock, P. J., Franzen, T. M. O., et al. 2019b, PASA, 36, e047
  * Kang & Koo (2007) Kang, J.-h. & Koo, B.-C. 2007, ApJS, 173, 85
  * Kang et al. (2012) Kang, J.-h., Koo, B.-C., & Salter, C. 2012, AJ, 143, 75
  * Khabibullin et al. (2023) Khabibullin, I. I., Churazov, E. M., Bykov, A. M., Chugai, N. N., & Sunyaev, R. A. 2023, MNRAS, 521, 5536
  * Koo & Kang (2004) Koo, B.-C. & Kang, J.-h. 2004, MNRAS, 349, 983
  * Koo et al. (2006) Koo, B.-C., Kang, J.-h., & Salter, C. J. 2006, ApJ, 643, L49
  * Koo et al. (2020) Koo, B.-C., Kim, C.-G., Park, S., & Ostriker, E. C. 2020, ApJ, 905, 35
  * Lu et al. (2021) Lu, C.-Y., Yan, J.-W., Wen, L., & Fang, J. 2021, Research in Astronomy and Astrophysics, 21, 033
  * McKee & Ostriker (1977) McKee, C. F. & Ostriker, J. P. 1977, ApJ, 218, 148
  * Pavlović et al. (2018) Pavlović, M. Z., Urošević, D., Arbutina, B., et al. 2018, ApJ, 852, 84
  * Reich & Sun (2019) Reich, W. & Sun, X.-H. 2019, Research in Astronomy and Astrophysics, 19, 045
  * Reid et al. (2016) Reid, M. J., Dame, T. M., Menten, K. M., & Brunthaler, A. 2016, ApJ, 823, 77
  * Reid et al. (2019) Reid, M. J., Menten, K. M., Brunthaler, A., et al. 2019, ApJ, 885, 131
  * Ripple et al. (2013) Ripple, F., Heyer, M. H., Gutermuth, R., Snell, R. L., & Brunt, C. M. 2013, MNRAS, 431, 1296
  * Shan et al. (2012) Shan, W. L., Yang, J., Shi, S. C., et al. 2012, IEEE Transactions on Terahertz Science and Technology, 2, 593
  * Sofue (2021) Sofue, Y. 2021, Galaxies, 9, 13
  * Sofue et al. (2021) Sofue, Y., Kohno, M., & Umemoto, T. 2021, ApJS, 253, 17
  * Stil et al. (2006) Stil, J. M., Taylor, A. R., Dickey, J. M., et al. 2006, AJ, 132, 1158
  * Su et al. (2017) Su, Y., Zhou, X., Yang, J., et al. 2017, ApJ, 845, 48
  * Tammann et al. (1994) Tammann, G. A., Loeffler, W., & Schroeder, A. 1994, ApJS, 92, 487
  * Toledo-Roy et al. (2014) Toledo-Roy, J. C., Velázquez, P. F., Esquivel, A., & Giacani, E. 2014, MNRAS, 437, 898
  * Wang et al. (2020) Wang, Y., Beuther, H., Rugel, M. R., et al. 2020, A&A, 634, A83
  * Wang et al. (2018) Wang, Y., Bihr, S., Rugel, M., et al. 2018, A&A, 619, A124
  * Wayth et al. (2015) Wayth, R. B., Lenc, E., Bell, M. E., et al. 2015, PASA, 32, e025
  * Wenger et al. (2000) Wenger, M., Ochsenbein, F., Egret, D., et al. 2000, A&AS, 143, 9
  * Xiao & Zhu (2014) Xiao, L. & Zhu, M. 2014, MNRAS, 438, 1081
  * Zhou et al. (2023) Zhou, X., Su, Y., Yang, J., et al. 2023, ApJS, 268, 61
  * Zhou et al. (2016) Zhou, X., Yang, J., Fang, M., et al. 2016, ApJ, 833, 4

##  Appendix A The schematic diagram

Figure 5: Schematic view of formation of expanding shell of G16.11$-$0.51. The
red star denotes the location of the initial SN explosion. The approximate
location of SNR G16.0$-$0.5 is also marked with a red circle. The arrows
indicate the velocities in different parts of the shell. According to the
distribution of the velocity component in the line of sight of different
sections of the shell, a similar structure will be present in the position-
velocity map as shown in Figure 1.

Figure 6: 12CO (J=1–0) intensity maps integrated over each 3 km s-1. The black
dashed circles are the same as in Figure 1. Central velocities are indicated
in each panel. The minimum value of each map is $3\sigma$.

We present a schematic diagram in Figure 5 to illustrate the formation of the
expanding shell of G16.11$-$0.51. As shown in the schematic diagram, when an
SNR breaks out of its parent MC, its shock wave would sweep up different
amounts of molecular gas in different directions with different decelerations,
resulting in an offset bubble structure. Such a structure will be further
enhanced by the radial density gradient of the MC. According to the
distribution of the velocity component in the line of sight of different parts
of the shell, a similar structure will be present in the position-velocity map
as shown in Figure 1. The blowout morphology has been observed in many SNRs,
e.g., 3C 400.2, RCW 103, and G352.7$-$0.1, which can be explained by the model
in which the supernova explosion producing the remnant occurs in the medium
with a density gradient or inside and near the border of a dense cloud
(Ambrocio-Cruz et al. 2006; Toledo-Roy et al. 2014; Lu et al. 2021). Figure 6
shows the distribution of 12CO (J=1–0) emission in the velocity range of $+$34
to $+$82 km s-1 with an interval of 3 km s-1.