11institutetext: University of California, Irvine, CA, USA 22institutetext:
University of Pennsylvania, PA, USA 33institutetext: The University of British
Columbia, BC, Canada

# On the Fairness of Swarm Learning in Skin Lesion Classification

Di Fan 11 Yifan Wu 22 Xiaoxiao Li Corresponding Author: Xiaoxiao Li,
<EMAIL_ADDRESS>
###### Abstract

Fairness is essential for trustworthy Artificial Intelligence (AI) in
healthcare. However, the existing AI model may be biased in its decision
marking. The bias induced by data itself, such as collecting data in subgroups
only, can be mitigated by including more diversified data. Distributed and
collaborative learning is an approach to involve training models in massive,
heterogeneous, and distributed data sources, also known as nodes. In this
work, we target on examining the fairness issue in Swarm Learning (SL), a
recent edge-computing based decentralized machine learning approach, which is
designed for heterogeneous illnesses detection in precision medicine. SL has
achieved high performance in clinical applications, but no attempt has been
made to evaluate if SL can improve fairness. To address the problem, we
present an empirical study by comparing the fairness among single (node)
training, SL, centralized training. Specifically, we evaluate on large public
available skin lesion dataset, which contains samples from various subgroups.
The experiments demonstrate that SL does not exacerbate the fairness problem
compared to centralized training and improves both performance and fairness
compared to single training. However, there still exists biases in SL model
and the implementation of SL is more complex than the alternative two
strategies.

##  1 Introduction

The success of deep learning can be partially attributed to data-driven
methodologies that automatically recognize patterns in large amounts of data.
Machine learning can theoretically be performed locally if adequate data and
computer infrastructure are available. However, data collection is costly and
data sharing, such as cloud computing, is strictly regulated due to the
privacy concerns in the healthcare domain.

Distributed and collaborative learning has gained popularity as a viable
solution to federate disparate data and computational resources to provide a
unifying analysis platform. Federated learning (FL) [15, 23] is a trending
type of learning scheme that avoids centralizing data in model training. Local
data owners (also known as nodes or clients) can train the private model
locally before sending the model weights or gradients to the central server
via FL. The central server then aggregates the shared model parameters to
create a new global model, which it then broadcasts to each local client.
Recently, a new format of collaborative learning has been proposed, called
Swarm Learning (SL) [21]. Unlike FL, SL replaces the central server to
coordinate model updates and parameters communication with the Swarm network,
which secures data privacy through blockchain technique [17] (See Figure
1(c)). Additionally, in SL, transactions can only be carried out by the nodes
having been pre-authorized and new nodes can be enrolled dynamically through
blockchain smart contract. Model parameters are exchanged in the Swarm network
by encryption and aggregated to update a new model at each round. We will
illustrate more on the SL system in Section 3.

The statistical heterogeneity in collaborative learning and the local model
multi-step training and then aggregation-based distributed learning paradigm
leaves the deep learning model vulnerable to biases in the model, as the model
may only attach importance to the dominant subgroup(s), such as, user groups
with sensitive (protected) attributes (i.e. age or sex) that are over-
represented in the deep learning model [7, 4]. Such fairness problems have
been recently discussed in the federated learning framework [12]. However, as
SL is a fresh new concept proposed in recent days, no attempt has been made to
evaluate the fairness or model bias issues in SL. Knowing the performance gap
in different subgroups can guide researchers and practitioners to seek methods
for better performance and fairness trade-off. This motivates us to conduct
evaluations on the fairness in medical imaging tasks with different models.
Specifically, we aim to compare model bias in centralized learning (pooling
data to a data center), SL, and training model with a subgroup of data only
(Figure 1). Without additional bias mitigation methods, we explore whether the
fairness behavior of SL leans to centralized learning or training on the
subgroup. To this end, we conduct experiments on the skin lesion
classification task, as skin lesions occur in diverse populations and have a
severe class imbalance.

Our highlights are summarized in threefold:

  1. 1.

We implement SL for a new medical application — skin lesion classification.

  2. 2.

To the best of our knowledge, we are the first to investigate fairness on SL.

  3. 3.

We observe that SL is robust to heterogeneous demographic-specific data
distributions on our task, and it does not degrade the performance and
fairness of the model compared to classical centralized training.

Figure 1: Different training strategies to be evaluated in this study.

##  2 Related Works

###  2.1 Collaborative learning and their application on healthcare

Collaborative learning is a decentralized/distributed system that uses the
principle of remote execution, which entails distributing duplication of a
machine learning algorithm to the sites where the data is stored. These sites
or devices are also called nodes. The system then runs the training method
locally before sending the results to a central repository to update the
algorithm [9]. In the development of AI algorithms in the healthcare domain,
data is one of the most important factors. However, there are two long-
standing data-related challenges: 1) There are no standardized patient
records; 2) Patient data is subject to stringent rules and protection
standards. As a result, the concept of federated learning has recently gained
a lot of attention. In the field of medical imaging, federated learning’s
applicability and benefits have been explored, such as on the problem of
whole-brain segmentation and brain tumor segmentation [16], as well as finding
disease-related biomarkers on fMRI [14].

###  2.2 Security and privacy of federated learning

In federated learning, security and privacy are critical, for example, in the
communication process between nodes. The existing methods include 1.
Differential privacy [6], which is a method of preserving a dataset’s global
statistical distribution while eliminating individually identifiable data; 2.
Homomorphic encryption, which is a type of encryption that allows unencrypted
data to be processed as if it were unencrypted data; 3. Secure multi-party
computation [24], meaning the process is done on encrypted data shares that
are distributed among them in such a way that no single party can get all of
the data on their own.

###  2.3 Fairness

Data inherently is biased. Deep learning algorithms have a tendency to mimic,
if not increase, data bias. For example, existing chest X-ray classifiers were
discovered to have differences in performance between subgroups defined by
distinct phenotypes [18]. Using gender-balanced datasets, a recent study found
statistically significant disparities in performance on medical imaging-based
diagnosis [11]. Darker skinned patients may be underrepresented [10, 13] in
existing dermatology datasets [20, 3], similar to forms of prejudice
associated to face recognition [22]. The effects of decision-making that is
(partially) based on the values of biased qualities can be irrevocable or even
lethal, especially in medical applications.

##  3 Problem Setting and Methods

###  3.1 Problem Setting

In this paper, we focus on skin lesion classification. Skin cancer is a server
public healthcare problem. In the United States, over 5 million newly
diagnosed cases every year [8]. Early detection of skin lesions will decrease
the death rate and reduce medical costs. Dermoscopy is one of the most widely
used skin imaging techniques to distinguish the lesion spots on skin [2]. Our
classification task is based on dermoscopic images.

Specifically, we cluster the dermoscopic images into subgroups based on
patients’ age and sex, which mimic the demographic-specific data distributions
in different data silos. We consider data is distributed and assume each
institution only contains a certain subgroup of dermoscopic images. We
evaluate the model performance and fairness on unseen testing data set. To
investigate model performance and bias in SL, we conduct experiments on skin
lesion classification tasks in SL, and compare it to centralized training
(pulling data together) and signal node training.

Next, we will introduce the different training strategies.

###  3.2 Swarm Learning

Swarm Learning(SL) is a decentralized, privacy-preserving Machine Learning
framework. It utilizes a dedicated server to train the Machine Learning models
with distributed data sources in a blockchain-safe strategy. As shown in
Figure 1(c), a set of local nodes process their local training data
respectively without sharing with each other to obtain an ML model
collaboratively. Then, via a Swam network rather than a centralized server,
the ML parameters or weights are shared by each individual node, and finally,
a merging model is formed.

At the beginning of SL, each local node should enroll or register with a
blockchain Swarm smart contract. Then, this one-time process could enable each
node to record information from the contract, like obtaining the model and
performing local training as soon as it meets the synchronization condition.
Next, the Swarm application programming interface (API) makes parameters
exchanging among nodes and merges to update a new model before a new round of
training. Finally, after the final round of parameter sharing and merging, the
framework will check the stop criterion and halt if the criterion is reached
otherwise back to local training.

With the blockchain techniques, SL has the following characteristics and
advantages: (1) Saving large medical data locally; (2) No need to exchange
original data, thereby reducing data traffic; (3) Providing high-level data
security and protecting model from attacks; (4) No need for a secure central
network; and (5) Allowing all members to merge parameters with equal rights.

###  3.3 Local and Centralized Training

Now, we describe two conventional training strategies to be compared with SL.

Single (Node) Training assumes no inter-institutional model and data sharing.
The institution uses the local data to train deep learning models.

Centralized Training collects the distributed data in the central database or
cloud server and trains deep learning models following the classical approach.

###  3.4 Fairness Definition and Metrics

In this section, we introduce the notation of fairness in classification model
and the fairness evaluation metrics to be used in this study.

######  Definition 1 (Fair classifier)

Let denote binarized sensitive attribute $z\in\\{0,1\\}$ that induces bias or
unfairness. We define a classifier $f$ with respect to data distribution $P$
on $\\{(X,Z,Y):\mathbb{R}^{d}\times\mathcal{Z}\times\mathcal{C}\\}$, where
$\mathcal{C}=\\{1,2,\dots,C\\}$ and $C$ is the number of classes. The
classifier $f$ is fair if

$\displaystyle P(f(x,z)\in\mathcal{C}|z=0)=P(f(x,z)\in\mathcal{C}|z=1)$

Namely, if the distance between the two output distributions of $z=0$ and
$z=1$ is zero, then $f$ is claimed as a fair classifier.

Before introducing the fairness metrics, we first give the formula of true
positive rate (TPR) and false positive rate (FPR). Given $x_{i},y_{i},z_{i}$
as input, label, and bias indicator for the $i$th sample in the dataset
$\mathcal{D}=\\{X,Y,Z\\}$. We denote the data samples in class
$c\in\mathcal{C}$ as $\mathcal{D}_{c}$. Let $\hat{y}_{i}$ represent the
predicted label of sample $i$.

The widely used bias quantification metrics are Statistical Parity Difference
(SPD), and Equal opportunity difference (EOD) [1, 5, 7]:

#### SPD

measures the difference in the probability of positive outcome between the
privileged and under-privileged groups:

$\displaystyle
SPD=\frac{1}{C}\sum_{c}^{C}\left(P_{(x_{i},y_{i},z_{i})\in\mathcal{D}_{c}}(\hat{y}_{i}=y_{i}|z_{i}=0)-P_{(x_{i},y_{i},z_{i})\in\mathcal{D}_{c}}(\hat{y}_{i}=y_{i}|z_{i}=1)\right).$

#### EOD

measures the difference in TPR for the privileged and under-privileged groups.
TPR for a given class $c\in\mathcal{C}$ and subgroup $z$ are defined as
$TPR_{z}^{c}=P_{(x_{i},y_{i},z_{i})\in\mathcal{D}_{c}}(\hat{y}_{i}=y_{i}|z_{i}=z)$.

$\displaystyle
EOD=\frac{1}{C}\sum_{c}^{C}\left(TPR_{z=0}^{c}-TPR_{z=1}^{c}\right).$

##  4 Experiment and Results

###  4.1 Dataset

Our data comes from the public skin lesion analysis dataset, Skin ISIC 2018
[20, 3]. However, we will divide the data into distributed portions to mimic
multi-institutional data collection. The whole ISIC 2018 dataset consists of
327 actinic keratosis (AKIEC), 514 basal cell carcinoma (BCC), 115
dermatofibroma (DF), 1113 melanoma (MEL), 6705 nevus (NV), 1099 pigmented
benign keratosis (BKL), 142 vascular lesions (VASC) smaples and in total 10015
RGB images. We resize the images to $224\times 224$. We consider sensitive
attributes of ISIC dataset are age ( $\geq$60 and $<$60), and sex (male and
female).

To evaluate our algorithms, we randomly split the dataset into a training set
(80%) and a test set (20%). As the data is unbalanced across lesion classes,
we augment $15\times,10\times,5\times,50\times,0\times,40\times,5\times$ for
AKIEC, BCC, DF, MEL, NV, BKL, and VASC, respectively. We divide the training
data into four subgroups – Male and age $\geq$ 60; Male and age $<$ 60; Female
and age $\geq$ 60; Female and age $<$ 60\. Hence we model four institutions,
and each owns the data from a subgroup. After augmentation, the population
proportion of the four subgroups are shown in Figure 2(a). For testing
dataset, to evaluate the model more efficiently, we make the data nearly
satisfy the condition of age $<$ 60: age $>$ 60 = 1:1 and male: female = 1:1
without augmentation over the four subgroups are shown in Figure 2(b).
Finally, figure 3 presents the heterogeneous sample distributions of the seven
lesion types in the four subgroups (here, the four institutions).

(a)

(b)

Figure 2: Training data distribution (left), Test data distribution (right).
The subgroup type is denoted next to the pie chart, where ‘M’ means ‘male’,
‘F’ means ‘female’, and ages are separated by $<60$ and $\geq 60$.

(a) Sample class distribution of M, $<$60

(b) Sample class distribution of M, $>=$60

(c) Sample class distribution of F, $<$60

(d) Sample class distribution of F, $>=$60

Figure 3: Label distributions in the four subgroups (institutions).

###  4.2  Implementation details

We use VGG19[19] with cross-entropy loss to classify skin lesions among the 7
classes (AKIEC, BCC, DF, MEL, NV, BKL, and VASC). Batch size is set to 32. We
use Adam optimizer, learning rate (lr) for classification training is 1e3 for
first 10 epochs, and 1e4 after 10 epochs.

For SL, we use a cluster of V100 GPUs to achieve our SL framework implemented
on Tensorflow 2.2. Eight individual IP addresses correspond to one License
Server, one Spire Server, two Swarm Network nodes, and four Swarm Learning
Nodes. Here, four Swarm learning Nodes represent four hospitals. Each Swarm
learning node has four kinds of data, respectively. All Swarm Network and
Swarm Learning nodes connect to the License Server, apls, running on host
172.17.0.3. All Swarm Network and Swarm Learning nodes connect to a single
SPIRE Server running on host 172.17.0.2. Two Swarm Network nodes utilized are
sn-1 and sn-2. sn-1 is the Sentinel Node running on host 172.17.0.4, and sn-2
runs on host 172.17.0.5. Each Swarm Network node has two Swarm Learning nodes
connected to it — Swarm Learning nodes sl-1 and sl-2 connect to sn- 1 while
sl-3 and sl-4 connect to sn-2. sl-1 runs on host 172.17.0.6, sl-2 runs on host
172.17.0.7, sl-3 runs on host 172.17.0.8, sl-4 runs on host 1172.17.0.9.

For single training, model is trained on one V100 GPU using only one subgroup
of data (M, age$<60$) representing a hospital.

For centralized training, we use the whole training dataset by pooling the
four hospitals in SL together and train the data with one V100 GPU.

Table 1: Model performance. Higher scores indicate better performance.

| Sex | Precision | Recall | F1-score | Age | Precision | Recall | F1-score
---|---|---|---|---|---|---|---|---
Centralized | Male | 0.774 | 0.680 | 0.707 | $\geq 60$ | 0.812 | 0.752 | 0.776
Female | 0.899 | 0.893 | 0.892 | $<60$ | 0.839 | 0.900 | 0.834
SL | Male | 0.691 | 0.658 | 0.656 | $\geq 60$ | 0.760 | 0.698 | 0.712
Female | 0.758 | 0.731 | 0.737 | $<60$ | 0.638 | 0.762 | 0.660
Single | Male | 0.579 | 0.460 | 0.498 | $\geq 60$ | 0.578 | 0.428 | 0.471
Female | 0.786 | 0.547 | 0.617 | $<60$ | 0.761 | 0.723 | 0.716

Table 2: Fairness scores. Lower SPD and EOD indicate less bias.

| Sex | Age
---|---|---
| Centralized |  SL | Single | Centralized |  SL | Single
SPD | 0.125 |  0.067 | 0.207 | 0.027 |  0.122 | 0.183
EOD | 0.213 |  0.073 | 0.087 | 0.148 |  0.064 | 0.295

###  4.3 Biases in Models Trained with Different Strategies

We evaluate our framework on skin lesion detection against sensitive
attribute: age and sex. The testing dataset contains samples from all the
above four subgroups. We present the bias that existed in the centralized
training, single training, and SL by showing the performance metrics
(precision, recall, F1-score). The results of centralized training are listed
in the ‘Centralized’ row of Table 1. Row ‘SL’ and row ‘Single’ list the
results for SL and single training respectively. The averaged precision,
recall, F1-score are reported to evaluate the testing classification
performance on the subgroups with different sensitive attributes. Centralized
training achieved the best results. The result is within our expectation as
the distributed optimization method in SL may hurt model performance,
especially on heterogeneous data. Nevertheless, SL still outperforms single
training and generalizes much better to testing data that contains unseen
subgroups, on the average performances of both subgroups.

We notice the performance gaps between the subgroups exist for all the
investigated training strategies, which reveals the biases. We denote the
group with better classification performance (such as female, $\rm age<60$) as
privilege group $z=0$ and the opposite group (such as male, $\rm age\geq 60$)
as under-privileged group $z=1$. We quantitatively measure the fairness scores
(SPD, EOD) on the testing set for our models and the vanilla model, as shown
in Table 2. Lower fairness scores indicate the model is less biased. We
observe that SL achieves best results for fairness scores on ‘sex’ partition
and EOD for ‘age’ partition, whereas the SPD for ‘age’ partition is less than
centralized. In this case, SL achieves comparable performance and the same
level of fairness in the heterogeneous demographic-specific data setting as
centralized training. Also, centralized has a higher EOD for ‘sex’ partition
than single and this could be in an explanation that single training has a
poor performance for both ‘male’ and ‘female’ recall value in Table 1 and
causes the difference to be small thus makes EOD lower compared to
centralized.

##  5 Discussion and Conclusion

In this work, we make the first attempt to evaluate the model fairness in SL
for medical applications. Specifically, we compare model performance and
fairness with various metrics on single, centralized, and SL. Based on our
experimental observation, SL could achieve better performance than training on
a single institution, and the SL model does not amplify biases. However, we
cannot conclude that SL gains the best performance-fairness trade-off for
arbitrary clinical tasks. We also want to point out the high implementation
complexity of SL framework due to the intricate blockchain network
configurations. In summary, our work brings the interesting collaborative
learning formula to skin cancer classification and designs experiments to
evaluate its fairness medical applications and highlight. Our future work
includes improving model performance for SL, investigating the mechanism of SL
in handling model fairness, and designing bias mitigation methods for SL.

## References

  * [1] Bellamy, R.K., Dey, K., Hind, M., Hoffman, S.C., Houde, S., Kannan, K., Lohia, P., Martino, J., Mehta, S., Mojsilovic, A., et al.: Ai fairness 360: An extensible toolkit for detecting, understanding, and mitigating unwanted algorithmic bias. arXiv preprint arXiv:1810.01943 (2018)
  * [2] Binder, M., Schwarz, M., Winkler, A., Steiner, A., Kaider, A., Wolff, K., Pehamberger, H.: Epiluminescence microscopy: a useful tool for the diagnosis of pigmented skin lesions for formally trained dermatologists. Archives of dermatology 131(3), 286–291 (1995)
  * [3] Codella, N.C., Gutman, D., Celebi, M.E., Helba, B., Marchetti, M.A., Dusza, S.W., Kalloo, A., Liopyris, K., Mishra, N., Kittler, H., et al.: Skin lesion analysis toward melanoma detection: A challenge at the 2017 international symposium on biomedical imaging (isbi), hosted by the international skin imaging collaboration (isic). In: 2018 IEEE 15th International Symposium on Biomedical Imaging (ISBI 2018). pp. 168–172. IEEE (2018)
  * [4] Du, M., Yang, F., Zou, N., Hu, X.: Fairness in deep learning: A computational perspective. IEEE Intelligent Systems pp. 1–1 (2020)
  * [5] Dwork, C., Hardt, M., Pitassi, T., Reingold, O., Zemel, R.: Fairness through awareness. In: Proceedings of the 3rd innovations in theoretical computer science conference. pp. 214–226 (2012)
  * [6] Dwork, C., Roth, A., et al.: The algorithmic foundations of differential privacy. Found. Trends Theor. Comput. Sci. 9(3-4), 211–407 (2014)
  * [7] Hardt, M., Price, E., Srebro, N.: Equality of opportunity in supervised learning. Advances in neural information processing systems 29, 3315–3323 (2016)
  * [8] Jerant, A.F., Johnson, J.T., Sheridan, C.D., Caffrey, T.J.: Early detection and treatment of skin cancer. American family physician 62(2), 357–368 (2000)
  * [9] Kaissis, G.A., Makowski, M.R., Rückert, D., Braren, R.F.: Secure, privacy-preserving and federated machine learning in medical imaging. Nature Machine Intelligence 2(6), 305–311 (2020)
  * [10] Kinyanjui, N.M., Odonga, T., Cintas, C., Codella, N.C.F., Panda, R., Sattigeri, P., Varshney, K.R.: Fairness of classifiers across skin tones in dermatology. In: Martel, A.L., Abolmaesumi, P., Stoyanov, D., Mateus, D., Zuluaga, M.A., Zhou, S.K., Racoceanu, D., Joskowicz, L. (eds.) Medical Image Computing and Computer Assisted Intervention – MICCAI 2020. pp. 320–329. Springer International Publishing, Cham (2020)
  * [11] Larrazabal, A.J., Nieto, N., Peterson, V., Milone, D.H., Ferrante, E.: Gender imbalance in medical imaging datasets produces biased classifiers for computer-aided diagnosis. Proceedings of the National Academy of Sciences 117(23), 12592–12594 (2020)
  * [12] Li, T., Sanjabi, M., Beirami, A., Smith, V.: Fair resource allocation in federated learning. arXiv preprint arXiv:1905.10497 (2019)
  * [13] Li, X., Cui, Z., Wu, Y., Gu, L., Harada, T.: Estimating and improving fairness with adversarial learning. arXiv preprint arXiv:2103.04243 (2021)
  * [14] Li, X., Gu, Y., Dvornek, N., Staib, L.H., Ventola, P., Duncan, J.S.: Multi-site fmri analysis using privacy-preserving federated learning and domain adaptation: Abide results. Medical Image Analysis 65, 101765 (2020)
  * [15] McMahan, B., Moore, E., Ramage, D., Hampson, S., y Arcas, B.A.: Communication-efficient learning of deep networks from decentralized data. In: Proceedings of Artificial Intelligence and Statistics (AISTATS). pp. 1273–1282. PMLR (2017)
  * [16] Rieke, N., Hancox, J., Li, W., Milletari, F., Roth, H.R., Albarqouni, S., Bakas, S., Galtier, M.N., Landman, B.A., Maier-Hein, K., et al.: The future of digital health with federated learning. NPJ digital medicine 3(1), 1–7 (2020)
  * [17] Saito, K., Yamada, H.: What’s so different about blockchain?—blockchain is a probabilistic state machine. In: 2016 IEEE 36th International Conference on Distributed Computing Systems Workshops (ICDCSW). pp. 168–175. IEEE (2016)
  * [18] Seyyed-Kalantari, L., Liu, G., McDermott, M., Chen, I.Y., Ghassemi, M.: Chexclusion: Fairness gaps in deep chest x-ray classifiers (2020)
  * [19] Simonyan, K., Zisserman, A.: Very deep convolutional networks for large-scale image recognition. arXiv preprint arXiv:1409.1556 (2014)
  * [20] Tschandl, P., Rosendahl, C., Kittler, H.: The ham10000 dataset, a large collection of multi-source dermatoscopic images of common pigmented skin lesions. Scientific data 5, 180161 (2018)
  * [21] Warnat-Herresthal, S., Schultze, H., Shastry, K.L., Manamohan, S., Mukherjee, S., Garg, V., Sarveswara, R., Händler, K., Pickkers, P., Aziz, N.A., et al.: Swarm learning for decentralized and confidential clinical machine learning. Nature 594(7862), 265–270 (2021)
  * [22] Wu, Y., Yang, F., Xu, Y., Ling, H.: Privacy-protective-gan for privacy preserving face de-identification. Journal of Computer Science and Technology (2019)
  * [23] Yang, Q., Liu, Y., Chen, T., Tong, Y.: Federated machine learning: Concept and applications. ACM Transactions on Intelligent Systems and Technology (TIST) 10(2), 12 (2019)
  * [24] Zhao, C., Zhao, S., Zhao, M., Chen, Z., Gao, C.Z., Li, H., Tan, Y.a.: Secure multi-party computation: theory, practice and applications. Information Sciences 476, 357–372 (2019)


# Decidability and $k$-Regular Sequences

Daniel Krenn  Fachbereich Mathematik, Paris Lodron University of Salzburg,
Hellbrunnerstraße 34, 5020 Salzburg, Austria<EMAIL_ADDRESS>and  Jeffrey
Shallit  School of Computer Science, University of Waterloo, Waterloo, ON N2L
3G1, Canada<EMAIL_ADDRESS>
###### Abstract.

In this paper we consider a number of natural decision problems involving
$k$-regular sequences. Specifically, they arise from considering

  * •

lower and upper bounds on growth rate; in particular boundedness,

  * •

images,

  * •

regularity (recognizability by a deterministic finite automaton) of preimages,
and

  * •

factors, such as squares and palindromes,

of such sequences. We show that these decision problems are undecidable.

Daniel Krenn is supported by the Austrian Science Fund (FWF): P 28466-N35.

This work was partially carried out by the authors during the workshop
“Numeration and Substitution 2019” at the Erwin Schrödinger International
Institute for Mathematics and Physics (ESI), and they were supported by it.

This paper is published under a Creative Commons Attribution 4.0 International
License.  (http://creativecommons.org/licenses/by/4.0/) Copyright held by the
authors.

##  1\. Introduction

A sequence $(a(n))_{n\geq 0}$ over a finite alphabet is said to be
$k$-automatic, for $k\geq 2$ an integer, if its $k$-kernel

$K_{k}({(a(n))_{n\geq 0}})=\\{(a(k^{e}n+i))_{n\geq 0}\colon\mathopen{}e\geq
0,\ 0\leq i<k^{e}\\}$

is of finite cardinality. There are many different equivalent definitions of
this class of sequences [1]. It is well-known that many questions about these
sequences, such as the growth rate of $\sum_{0\leq n<N}a(n)$, are decidable
[5].

The so-called $k$-regular sequences form a natural generalization of the
automatic sequences. These are sequences $(a(n))_{n\geq 0}$ where the kernel
$K_{k}((a(n))_{n\geq 0})$ is contained in a finitely generated module. Unlike
the case of $k$-automatic sequences, it is known that some decision problems
involving $k$-regular sequences are recursively unsolvable [2].

In this paper we examine a number of natural decision problems involving
$k$-regular sequences, and show that they are undecidable.

Some general results have recently been obtained by Honkala [11].

###  1.1. Recursively solvable decision problems

A decision problem is one with a yes/no answer. To say that a decision problem
is solvable means there exists an algorithm (or Turing machine) that
unerringly solves it on all inputs. Throughout this paper we use the terms
“recursively solvable”, “solvable”, and “decidable” interchangeably, and
similarly for the terms “recursively unsolvable”, “unsolvable”, and
“undecidable”.

###  1.2. Notation

We let $\mathbb{N}_{0}$ denote the nonnegative integers (natural numbers) and
$\mathbb{N}$ denote the positive integers.

For a word $z$ with symbols chosen from a finite set $D$, we let $\lvert
z\rvert$ denote its length and $\lvert z\rvert_{d}$ the number of occurrences
of the letter $d\in D$ in $z$.

For a fixed integer $k\geq 2$, we consider base-$k$ representations with the
usual digit set $D=\\{0,1,\dots,k-1\\}$. For a nonnegative integer $n$, we
write $(n)_{k}$ for the standard $k$-ary representation of $n$, most
significant digit first, and having no leading zeroes. The representation of
$0$ is the empty word. Note that $(n)_{k}$ is a word over $D$ and that
$\lvert(n)_{k}\rvert=\lfloor\log_{k}n\rfloor+1$ for $n>0$.

###  1.3. Hilbert’s tenth problem

Showing that a certain decision problem is recursively unsolvable is often
carried out by constructing a (Turing) reduction from another decision problem
already known to be recursively unsolvable. One such problem is Hilbert’s
tenth problem [6, 17]:

######  Theorem (Hilbert’s tenth problem; variant).

The decision problem

> “Given a multivariate polynomial $p$ with integer coefficients, do there
> exist natural numbers $x_{1}$, $x_{2}$, …, $x_{t}$ such that
> $p(x_{1},\dots,x_{t})=0$?”

is recursively unsolvable.

The analogous problem, where the $x_{i}$ need to be positive, is also
recursively unsolvable. We will reduce from this problem quite frequently,
namely in Theorems B, H, I, J and L.

###  1.4. Representation of $k$-regular sequences

A $k$-regular sequence $(f(n))_{n\geq 0}$ can be finitely represented in a
number of different ways, of which two are the most useful. First, a set of
identities in terms of sequences in the $k$-kernel, where each identity
represents a subsequence $(f(k^{e}n+i))_{n\geq 0}$, $0\leq i<k^{e}$, as a
linear combination of subsequences in the $k$-kernel, and a set of initial
values. Together it must be possible to compute $f(n)$ for all $n$ from this
set of identities and initial values.

Second, a linear representation for $(f(n))_{n\geq 0}$, which consists of a
$1\times r$ row vector $v$, an $r\times 1$ column vector $w$, and $k$ square
matrices $M_{0},M_{1},\ldots,M_{k-1}$ of dimension $r\times r$ such that

$f(n)=vM_{n_{0}}\cdots M_{n_{s-1}}w,$

for all $n$, where $(n)_{k}=n_{s-1}\cdots n_{0}$. Of course, the empty product
of matrices is the identity matrix. The rank of a linear representation is
$r$. See [1, Theorem 2.2]. A linear representation is minimal if it has
smallest possible rank for the corresponding sequence.

For example, consider the $2$-regular sequence $(s_{2}(n))_{n\geq 0}$, which
counts the number of $1$’s in the binary representation of $n$. Then it is
easy to see that

$\displaystyle s_{2}(0)$ $\displaystyle=0$ $\displaystyle s_{2}(2n)$
$\displaystyle=s_{2}(n)$ $\displaystyle s_{2}(4n+1)$
$\displaystyle=s_{2}(2n+1)$ $\displaystyle s_{2}(4n+3)$
$\displaystyle=-s_{2}(n)+2s_{2}(2n+1)$

is an example of the former representation, and

$v=[0\ \ 1];\quad M_{0}=\left[\begin{array}[]{cc}1&0\\\
0&1\end{array}\right];\quad M_{1}=\left[\begin{array}[]{cc}0&-1\\\
1&2\end{array}\right];\quad w=\left[\begin{array}[]{c}1\\\
0\end{array}\right]$

is an example of the latter.

From now on, when we say an algorithm is “given” a $k$-regular sequence as
input, we actually mean that the input is either one of these two
representations. Note that we can transform between these two representations
effectively, that is, with an algorithm [3].

Some of our theorems involve algebraic numbers. When we say we are “given” a
real algebraic number $\alpha$, we mean we are given the minimal polynomial
for $\alpha$, together with a rational interval that contains $\alpha$ and
none of its conjugates. As is well-known [9, 10], we can effectively carry out
arithmetic on algebraic numbers represented in this way.

###  1.5. Closure properties of $k$-regular sequences

In this section we recall some closure properties of $k$-regular sequences:
that is, which operations on sequences preserve the property of
$k$-regularity. For more details, see [1]. It is important to note that not
only do these operations preserve $k$-regularity; they also are effectively
$k$-regular. Let $\circ$ be some operation on sequences. By the operation
$\circ$ being “effectively $k$-regular”, we mean that there is an algorithm
that, given some representation of $k$-regular sequences ${\bf
a}=(a(n))_{n\geq 0}$ and ${\bf b}=(b(n))_{n\geq 0}$, computes a representation
for ${\bf a}\circ{\bf b}$.

######  Theorem.

The class of $k$-regular sequences is (effectively) closed under the following
operations:

  * (a)

sum, ${\bf a}+{\bf b}=(a(n)+b(n))_{n\geq 0}$;

  * (b)

product, ${\bf a}{\bf b}=(a(n)\,b(n))_{n\geq 0}$;

  * (c)

convolution, ${\bf a}\star{\bf b}=(\sum_{0\leq i\leq n}a(i)\,b(n-i))_{n\geq
0}$;

  * (d)

perfect shuffle, ${\bf a}\,\mathcyr{SH}\,{\bf b}={\bf c}=(c(n))_{n\geq 0}$,
where $c(2i)=a(i)$ and $c(2i+1)=b(i)$ for $i\geq 0$. The same is true for
$t$-way perfect shuffle, where we combine $t$ sequences analogously.

For proofs, see [1].

######  Remark 1.1.

Let $p$ be a multivariate polynomial with integer coefficients, and suppose
$d_{1},\ldots,d_{t}\in\\{0,\dots,k-1\\}$. Then

$(\mathop{{#1}{}}(\lvert z\rvert_{d_{1}},\lvert z\rvert_{d_{2}},\dots,\lvert
z\rvert_{d_{t}}))_{n\geq 0}$

with $z=(n)_{k}$, is (effectively) a $k$-regular sequence in
$n\in\mathbb{N}_{0}$ over $\mathbb{Z}$. This is true because the number of
occurrences $\lvert z\rvert_{d}$ of a digit $d$ in the standard $k$-ary
expansion $z=(n)_{k}$ is a $k$-regular sequence by [1, Theorem 6.1] and, as
above, $k$-regular sequences are closed under term-by-term addition and
multiplication.

##  2\. Equality of $k$-regular sequences

Before we consider growth of $k$-regular sequences in the next section, we
study decidability of whether one $k$-regular sequence equals another.

######  Theorem A.

Let $k\geq 2$ be an integer and $\mathbb{S}\supseteq\mathbb{N}_{0}$ be a set.
The decision problem

> “Given two $k$-regular sequences $(f(n))_{n\geq 0}$ and $(g(n))_{n\geq 0}$
> over $\mathbb{S}$, does $f(n)=g(n)$ for all $n$?”

is recursively solvable.

###### Proof.

The sequence $(f(n)-g(n))_{n\geq 0}$ is $k$-regular by the closure properties
of $k$-regular sequences. We compute a linear representation of this sequence
and apply a minimization algorithm (see Berstel and Reutenauer [3]) to it.
This results in a linear representation of minimum rank. This rank is $0$ iff
we started with the zero sequence, iff $f(n)=g(n)$ for all $n$. ∎

##  3\. Growth of $k$-regular sequences

We use the standard notation for asymptotic growth of sequences. Let
$(f(n))_{n\geq 0}$ and $(g(n))_{n\geq 0}$ be sequences. Then

  * •

$f(n)\in O(g(n))$ means that there exist $n_{0}$, $c>0$ such that $f(n)\leq
cg(n)$ for all $n\geq n_{0}$, and

  * •

$f(n)\in\Omega(g(n))$ means that there exist $n_{0}$, $c>0$ such that
$f(n)\geq cg(n)$ for all $n\geq n_{0}$.

For simplicity, we sometimes say that the sequence $(f(n))_{n\geq 0}$ is in
$O(g(n))$ or $\Omega(g(n))$.

In what follows, $\mathbb{A}$ denotes the set of real algebraic numbers.

###  3.1. Lower bounds

######  Theorem B.

Let $k\geq 2$ be an integer and $\mathbb{S}$ be a set with
$\mathbb{N}\subseteq\mathbb{S}\subseteq\mathbb{A}$. The decision problem

> “Given a $k$-regular sequence $(f(n))_{n\geq 0}$ over $\mathbb{S}$, is
> $f(n)$ in $\mathop{{#1}{}}(n)$?”

is recursively unsolvable.

###### Proof.

We reduce from Hilbert’s tenth problem. For a given multivariate polynomial
$p$ in $t$ variables over $\mathbb{Z}$, we choose $r\in\mathbb{N}$ such that
$K=k^{r}>t+1$, and we construct the sequence

$f(n)\coloneqq(n+1)\bigl{(}\mathop{{#1}{}}(\lvert z\rvert_{1},\lvert
z\rvert_{2},\dots,\lvert z\rvert_{t})\bigr{)}^{2}\,(\lvert z\rvert_{t+1}+1)$

with $z=(n)_{K}$. The sequence $(f(n))_{n\geq 0}$ is $K$-regular (see Remark
1.1) and therefore $k$-regular as well [1, Theorem 2.9].

The following claim shows that the above indeed provides a reduction.

######  Claim.

The sequence $(f(n))_{n\geq 0}$ is not in $\Omega(n)$ iff there exist
nonnegative integers $x_{1}$, $x_{2}$, …, $x_{t}$ such that
$p(x_{1},\dots,x_{t})=0$.

To see this, note that $f(n)=0$ iff at least one factor is zero, but the first
and third factors defining $f$ are never zero. Hence $f(n)=0$ iff
$\mathop{{#1}{}}(\lvert z\rvert_{1},\lvert z\rvert_{2},\dots,\lvert
z\rvert_{t})=0$. Moreover, note that if a zero of $(f(n))_{n\geq 0}$ occurs
once, it occurs infinitely often by the third factor of the definition for
$f$.

Thus, if $p(x_{1},\dots,x_{t})=0$ has at least one solution, then $f(n)=0$ for
infinitely many $n$, and hence $f(n)$ is not in $\Omega(n)$. Otherwise, if
$p(x_{1},\dots,x_{t})=0$ does not have any solution, then its absolute value
is at least $1$, so the value of $f(n)$ is at least $n+1$, and hence $f(n)$ is
in $\Omega(n)$. This completes the proof of the claim and consequently the
proof of Theorem B. ∎

Theorem B can be extended to other growth rates as well.

######  Corollary C.

Let $k\geq 2$ be an integer. Suppose $\sigma$ is a real number and $\ell$ is a
nonnegative integer, not both zero. Assume that $k^{\sigma}$ is an algebraic
number. Let $\mathbb{S}$ be a ring with
$\mathbb{Z}\subseteq\mathbb{S}\subseteq\mathbb{A}$ and containing
$k^{\sigma}$. The decision problem

> “Given a $k$-regular sequence $(f(n))_{n\geq 0}$ over $\mathbb{S}$, is
> $f(n)$ in $\mathop{{#1}{}}\big{(}n^{\sigma}(\log n)^{\ell}\big{)}$?”

is recursively unsolvable.

######  Remark 3.1.

For a real number $\sigma$ and a nonnegative integer $\ell$ we construct a
$k$-regular sequence $(h_{\sigma,\ell}(n))_{n\geq 0}$ with positive terms
(except for the first few terms, which may be $0$) satisfying

$h_{\sigma,\ell}(n)\in\mathop{{#1}{}}big(n^{\sigma}(\log n)^{\ell}big),$

as follows.

Set $H_{0}=\cdots=H_{k-1}=J_{\ell+1}(k^{\sigma})$, where
$J_{\ell+1}(k^{\sigma})$ is a Jordan block, of size $\ell+1$, corresponding to
the eigenvalue $k^{\sigma}$. We set

$h_{\sigma,\ell}(n)=e_{1}\mathop{{#1}{}}(n)\,e_{\ell+1},$

where $H(n)=H_{n_{0}}\cdots H_{n_{\ell-1}}$ for $(n)_{k}=n_{\ell-1}\cdots
n_{0}$, and the $e_{i}$ are the $i$th unit vectors. Therefore, the sequence
$h_{\sigma,\ell}(n)$ is $k$-regular, as it is defined by a linear
representation. Explicitly, we have

$h_{\sigma,\ell}(n)=\binom{s}{\ell}k^{(s-\ell)\sigma},$

where $s=\lfloor\log_{k}n\rfloor+1$. Thus, this sequence’s asymptotic behavior
is $\mathop{{#1}{}}big(n^{\sigma}(\log n)^{\ell}big)$. If $\ell=0$, then no
term is $0$. If $\ell\not=0$, then only terms with $n\leq k^{\ell-1}$ are $0$.

###### Proof of Corollary C.

The proof runs along the same lines as the proof of Theorem B. For a given
multivariate polynomial $p$ in $t$ variables over $\mathbb{Z}$, we instead
choose $r\in\mathbb{N}$ such that $K=k^{r}>t+1$, and we define

$f(n)\coloneqq(\mathop{{#1}{}}(n)+1)\bigl{(}\mathop{{#1}{}}(\lvert
z\rvert_{1},\lvert z\rvert_{2},\dots,\lvert z\rvert_{t})\bigr{)}^{2}\,(\lvert
z\rvert_{t+1}+1)$

with $z=(n)_{K}$ and $(\mathop{{#1}{}}(n))_{n\geq 0}$ of Remark 3.1. Note that
the factor $(\mathop{{#1}{}}(n)+1)$ of $f(n)$ is increasing and always
positive. ∎

###  3.2. Upper bounds

Let $(h(n))_{n\geq 0}$ be a sequence. We say that a sequence $(M(n))_{n\geq
0}$ of matrices with entries in a set $\mathbb{S}\subseteq\mathbb{A}$ is in
$O(h(n))$, formally written as usual as

$M(n)\in O(h(n)),$

if each sequence of a fixed entry (fixed row and column) of the matrices is in
$O(h(n))$. Rephrased, this means that the sequence of maximum norms of the
matrices lies in $O(h(n))$. By the equivalence of norms, this is also true for
any other norm. As in the one-dimensional case, we say that the sequence
$(M(n))_{n\geq 0}$ is _bounded_ if it lies in $O(1)$.

######  Remark 3.2.

Let $\sigma\in\mathbb{R}$ and $\ell\in\mathbb{N}_{0}$, and set
$h(n)=n^{\sigma}(\log n)^{\ell}$. If $m(n)=cn+\tau(n)$ for some constant
$c\neq 0$ and some sequence $\tau(n)\in o(n)$, then

$O(h(m))=O(h(n)).$

as $n\to\infty$. This follows from

$\displaystyle h(m(n))$
$\displaystyle=\bigl{(}cn+\tau(n)\bigr{)}^{\sigma}\bigl{(}\log(cn+\tau(n))\bigr{)}^{\ell}$
$\displaystyle=c^{\sigma}h(n)\Bigl{(}1+\frac{\tau(n)}{cn}\Bigr{)}^{\sigma}\biggl{(}1+\frac{\log\bigl{(}c+\frac{\tau(n)}{n}\bigr{)}}{\log
n}\biggr{)}^{\ell}$ $\displaystyle=c^{\sigma}h(n)\bigl{(}1+o(1)\bigr{)}.$

######  Theorem D.

Let $k\geq 2$ be an integer and $\mathbb{S}$ be a ring with
$\mathbb{Q}\subseteq\mathbb{S}\subseteq\mathbb{A}$. The decision problem

> “Given a $k$-regular sequence $(f(n))_{n\geq 0}$ over $\mathbb{S}$, is
> $f(n)$ bounded?”

is recursively unsolvable.

######  Remark 3.3.

The above problem is decidable for $k$-regular sequences that have a linear
representation with integer matrices; see the algorithm of Mandel and Simon
[16] for matrices with nonnegative entries and the algorithm of Jacob [13, 14,
15] for general integer matrices.

Given square matrices $F_{0}$, …, $F_{k-1}$ over a ring $\mathbb{S}$ all of
the same dimension, we set $F(n)=F_{n_{0}}\cdots F_{n_{s-1}}$ for
$(n)_{k}=n_{s-1}\cdots n_{0}$, and we call the sequence $(F(n))_{n\geq 0}$ the
matrix-valued linear representation sequence of the set
$\\{F_{0},\ldots,F_{k-1}\\}$. We use this notion in the following lemma, which
extends decidability from $k$-regular sequences to matrix-valued linear
representation sequences.

######  Lemma 3.4.

Let $\mathbb{S}$ be a ring. Let $P$ be a property of a sequence over
$\mathbb{S}$, i.e., $P$ for each sequence over $\mathbb{S}$ is either true or
false. Suppose we can extend property $P$ to sequences of matrices over
$\mathbb{S}$ in one of the following ways: Property $P$ holds for a sequence
of matrices iff $P$ holds for

  1. (1)

all sequences or

  2. (2)

any sequence

consisting of a fixed entry (fixed row and column).

If $P$ is recursively solvable for $k$-regular sequences over $\mathbb{S}$,
then $P$ is recursively solvable for matrix-valued linear representation
sequences of a set of $k$ square matrices over $\mathbb{S}$, all of the same
dimension.

We will use this lemma in the proof of Theorem D with the property $P$ being
the boundedness of a sequence and in the proof of Theorem F, where $P$ is true
iff a sequence does not have at least polynomial growth.

###### Proof of Lemma 3.4.

We show (1). Then (2) follows by using the negation of property $P$.

We define

$f_{i,j}(n)=e_{i}F(n)e_{j}$

with $F(n)$ as above and where $e_{i}$ is the $i$th unit vector. Therefore,
$(f_{i,j}(n))_{n\geq 0}$ is the sequence of entries in row $i$ and column $j$
of $(F(n))_{n\geq 0}$. All sequences $(f_{i,j}(n))_{n\geq 0}$ are $k$-regular,
as they are defined by a linear representation. Clearly all these sequences
$(f_{i,j}(n))_{n\geq 0}$ satisfy $P$ iff $F(n)$ satisfies $P$.

As the question of deciding property $P$ of a $k$-regular sequence is
recursively solvable, we can decide $P$ for each of the (finitely many)
distinct sequences $(f_{i,j}(n))_{n\geq 0}$ and therefore can decide $P$ for
$(F(n))_{n\geq 0}$. ∎

###### Proof of Theorem D.

We reduce from the question of boundedness of all products of matrices over
the rationals, which is not recursively solvable; see Blondel and Tsitsiklis
[4].

For a given set of matrices $\\{F_{0},\ldots,F_{k-1}\\}$, for all $n_{0}$, …,
$n_{s-1}\in\\{0,\dots,k-1\\}$, there is either a largest index
$j\in\\{1,\dots,s\\}$ with $n_{j-1}\neq 0$ or $j=0$ and such that
$F_{n_{0}}\cdots F_{n_{s-1}}=F_{n_{0}}\cdots F_{n_{j-1}}F_{0}^{s-j}$. Whenever
we are now deciding the boundedness of matrix products, we split as above and
consider the factors $F_{0}^{s-j}$ on the right-hand side and the remaining
product separately. As for the $F_{0}^{s-j}$, we can decide the boundedness of
powers of a single matrix from knowledge of its Jordan decomposition. And, as
for the remaining factors, the statement of the theorem follows by using the
reduction that is provided by Lemma 3.4 with property $P$ being the
boundedness of a sequence. ∎

######  Corollary E.

Let $k\geq 2$ be an integer, $\sigma$ a real number and $\ell$ a nonnegative
integer. Assume that $k^{\sigma}$ is an algebraic number. Let $\mathbb{S}$ be
a ring with $\mathbb{Q}\subseteq\mathbb{S}\subseteq\mathbb{A}$ and containing
$k^{\sigma}$. The decision problem

> “Given a $k$-regular sequence $(f(n))_{n\geq 0}$ over $\mathbb{S}$, is
> $f(n)$ in $O(n^{\sigma}(\log n)^{\ell})$?”

is recursively unsolvable.

###### Proof.

We reduce from the decision problem stated in Theorem D. For a $k$-regular
sequence $(g(n))_{n\geq 0}$, we construct
$\mathop{{#1}{}}(n)=\mathop{{#1}{}}(n)\mathop{{#1}{}}(n)$ with
$\mathop{{#1}{}}(n)$ as defined in Remark 3.1.

Then, the $k$-regular sequence $(f(n))_{n\geq 0}$ is in $O(n^{\sigma}(\log
n)^{\ell})$ iff $g(n)$ is in $O(1)$, i.e., bounded. Therefore deciding if a
$k$-regular sequence is in $O(n^{\sigma}(\log n)^{\ell})$ implies deciding the
boundedness of a $k$-regular sequence, which contradicts Theorem D. ∎

Let $\sigma\in\mathbb{R}$ and $\ell\in\mathbb{N}_{0}$. We say that a sequence
$(f(n))_{n\geq 0}$ has exact growth $n^{\sigma}(\log n)^{\ell}$ if

$f(n)\in O\big{(}n^{\sigma}(\log n)^{\ell}\big{)}$

but for all $\sigma^{\prime}\in\mathbb{R}$ and
$\ell^{\prime}\in\mathbb{N}_{0}$ with $(\sigma^{\prime},\ell^{\prime})$
lexicographically smaller than $(\sigma,\ell)$ we have

$f(n)\not\in O\big{(}n^{\sigma^{\prime}}(\log n)^{\ell^{\prime}}\big{)}.$

######  Proposition 3.5.

Let $(f(n))_{n\geq 0}$ be a $k$-regular sequence over a field
$\mathbb{S}\subseteq\mathbb{A}$ with matrices $F_{0},\dots,F_{k-1}$ of a
minimal linear representation, and set $F(n)=F_{n_{0}}\cdots F_{n_{s-1}}$ for
$(n)_{k}=n_{s-1}\cdots n_{0}$. Let $\sigma\in\mathbb{R}$ and
$\ell\in\mathbb{N}_{0}$, and set $h(n)=n^{\sigma}(\log n)^{\ell}$. Then

$f(n)\in O(h(n))$

if and only if

$F(n)\in O(h(n)).$

In particular, both $f(n)$ and $F(n)$ have exactly the same growth rate.

###### Proof.

Let $\lambda$ and $\gamma$ be the vectors of our minimal representation, i.e.,
$f(n)=\lambda\mathop{{#1}{}}(n)\gamma$ for all $n\in\mathbb{N}_{0}$. We start
with the easy direction: As $f(n)$ is a finite linear combination of the
entries in the matrix $F(n)$ and each of these entries is in $O(h(n))$, we
have $f(n)$ is in $O(h(n))$.

Conversely, suppose $F(n)$ is not in $O(h(n))$. As there are only finitely
many entries in each matrix $F(n)$, we can assume that one entry of $F(n)$ is
not in $O(h(n))$. Let $(g(n))_{n\geq 0}$ denote the sequence of this fixed
entry of the matrices.

Now, as our linear representation is minimal, there exist finite subsets $P$,
$Q\subseteq\mathbb{N}_{0}$ and coefficients $c_{p}$,
$d_{q}\in\mathbb{S}\setminus\\{0\\}$ for $p\in P$, $q\in Q$ such that

(3.1) $g(n)=\lambda\biggl{(}\,\sum_{p\in P}\sum_{q\in
Q}c_{p}d_{q}\mathop{{#1}{}}(p)\mathop{{#1}{}}(n)\mathop{{#1}{}}(q)\biggr{)}\gamma$

for all $n\in\mathbb{N}_{0}$; see [3, Corollary 2.3]. As $g(n)$ is not in
$O(h(n))$, one of the finitely many summands

$c_{p}d_{q}\cdot\lambda\mathop{{#1}{}}(p)\mathop{{#1}{}}(n)\mathop{{#1}{}}(q)\gamma$

of (3.1), where $p\in P$ and $q\in Q$, is not in $O(h(n))$. Dividing this
summand by $c_{p}d_{q}$ yields a subsequence of $(f(n))_{n\geq 0}$, namely

$\lambda\mathop{{#1}{}}(p)\mathop{{#1}{}}(n)\mathop{{#1}{}}(q)\gamma=\mathop{{#1}{}}(m(n))$

with
$m(n)=pk^{\lvert(n)_{k}\rvert+\lvert(q)_{k}\rvert}+nk^{\lvert(q)_{k}\rvert}+q$
for all $n\in\mathbb{N}_{0}$. As
$\lvert(n)_{k}\rvert=\lfloor\log_{k}n\rfloor+1$, we have $m(n)=cn+o(n)$ for
some constant $c$, and therefore, by Remark 3.2, we obtain that the
subsequence $\mathop{{#1}{}}(m(n))$ is not in $O(h(n))=O(h(m(n)))$. Thus the
sequence $(f(n))_{n\geq 0}$ itself is not in $O(h(n))$. ∎

###  3.3. Polynomial growth

The growth of a $k$-regular sequence is always at most polynomial. To be
precise, for a $k$-regular sequence $(f(n))_{n\geq 0}$ with values in
$\mathbb{A}$, there exists a real constant $\sigma\geq 0$ such that $f(n)\in
O(n^{\sigma})$; see [1, Theorem 2.10].

######  Theorem F.

Let $k\geq 2$ be an integer and $\mathbb{S}$ be a ring with
$\mathbb{Q}\subseteq\mathbb{S}\subseteq\mathbb{A}$. The decision problem

> “Given a $k$-regular sequence $(f(n))_{n\geq 0}$ over $\mathbb{S}$, does
> $f(n)$ have at least polynomial growth, i.e., does there exist $\sigma>0$
> such that $f(n)$ is not in $O(n^{\sigma})$?”

is recursively unsolvable.

Consequently, the decision problem in this theorem is whether a $k$-regular
sequence has polynomial growth or a smaller growth.

Recall that the joint spectral radius of a finite set $S$ of square matrices
is defined to be

$\rho(S)=\lim_{\ell\rightarrow\infty}\max\\{\lVert A_{1}\cdots
A_{\ell}\rVert^{1/\ell}\colon\mathopen{}A_{i}\in S\\},$

where $\lVert\,\cdot\,\rVert$ is any matrix norm.

######  Proposition 3.6.

Let $\mathbb{S}\subseteq\mathbb{A}$ be a ring, let $\rho\in\mathbb{R}$ be
positive, let $F_{0}$, …, $F_{k-1}$ be square matrices over $\mathbb{S}$ all
of the same dimension, and set $F(n)=F_{n_{0}}\cdots F_{n_{s-1}}$ for
$(n)_{k}=n_{s-1}\cdots n_{0}$. Then the following two statements are
equivalent:

  1. (1)

The joint spectral radius of $F_{0}$, …, $F_{k-1}$ is $\rho$.

  2. (2)

For all $\varepsilon>0$ we have $F(n)\in
O\big{(}n^{(\log_{k}\rho)+\varepsilon}\big{)}$ and $F(n)\not\in
O\big{(}n^{(\log_{k}\rho)-\varepsilon}\big{)}$ as $n\to\infty$, and we have
$F_{0}^{s}\in O\big{(}(k^{s})^{(\log_{k}\rho)+\varepsilon}\big{)}$ as
$s\to\infty$.

In particular, the joint spectral radius $\rho$ of $F_{0}$, …, $F_{k-1}$ is
bounded by some positive $\rho^{\prime}\in\mathbb{R}$, i.e.,
$\rho\leq\rho^{\prime}$, iff for all $\varepsilon>0$ we have $F(n)\in
O\big{(}n^{(\log_{k}\rho^{\prime})+\varepsilon}\big{)}$ as $n\to\infty$ and
$F_{0}^{s}\in O\big{(}(k^{s})^{(\log_{k}\rho^{\prime})+\varepsilon}\big{)}$ as
$s\to\infty$.

If $\mathbb{S}$ is a field and the matrices $F_{0}$, …, $F_{k-1}$ are of a
minimal representation of a $k$-regular sequence $(f(n))_{n\geq 0}$, then we
may replace $F(n)$ by $f(n)$ in the statements of this proposition.

######  Corollary G.

Let the assumptions of Proposition 3.6 hold and suppose $F_{0}$ is the zero-
matrix. Then the joint spectral radius $\rho$ of $F_{1}$, …, $F_{k-1}$ is
bounded by some positive $\rho^{\prime}\in\mathbb{R}$, i.e.,
$\rho\leq\rho^{\prime}$, iff for all $\varepsilon>0$ we have $F(n)\in
O\big{(}n^{(\log_{k}\rho^{\prime})+\varepsilon}\big{)}$ as $n\to\infty$.

We will use this corollary with $\rho^{\prime}=1$ in the proof of Theorem F to
connect polynomial growth with the joint spectral radius $\rho$.

###### Proof of Proposition 3.6.

In this proof, we suppose that $s$ and $n$ are related by
$s=\lfloor\log_{k}n\rfloor+1$. Then, by Remark 3.2, we have
$O(n^{\sigma})=O(k^{s\sigma})$ as $n\to\infty$ for any $\sigma$.

We have that for any fixed real $\sigma$,

$\lVert F(n)\rVert\in O(n^{\sigma})=O(k^{s\sigma})$

as $n\to\infty$ is equivalent to

(3.2) $\max_{k^{s-1}\leq n<k^{s}}\lVert F(n)\rVert\in O(k^{s\sigma})$

as $s\to\infty$, because $s$ is the same for all $n$ within the given range of
the argument of the maximum and the right-hand side $O(k^{s\sigma})$ only
depends on $s$ (and not on $n$).

We set

$\rho_{s}\coloneqq\max_{n_{0},\dots,n_{s-1}\in\\{0,\dots,k-1\\}}\lVert
F_{n_{0}}\cdots F_{n_{s-1}}\rVert^{1/s}$

Then the bound (3.2) together with $F_{0}^{s}\in O(k^{s\sigma})$ is equivalent
to

(3.3) $k^{s\log_{k}\rho_{s}}=\rho_{s}^{s}\in O(k^{s\sigma})$

as $s\to\infty$, because there is a constant $c>0$ (only depending on the used
norm) such that for all $n_{0}$, …, $n_{s-1}\in\\{0,\dots,k-1\\}$, there is
either a largest index $j\in\\{1,\dots,s\\}$ with $n_{j-1}\neq 0$ or $j=0$ and

$\lVert F_{n_{0}}\cdots F_{n_{s-1}}\rVert\leq c\lVert F_{n_{0}}\cdots
F_{n_{j-1}}\rVert\cdot\lVert F_{0}^{s-j}\rVert\in
O(k^{j\sigma})\,O(k^{(s-j)\sigma})=O(k^{s\sigma}).$

Consequently, the bound (3.3) is equivalent to the existence of an
$S\in\mathbb{N}_{0}$ such that for all $s\geq S$, the inequality
$\log_{k}\rho_{s}\leq\sigma$ holds. So much for our preliminary
considerations.

Now let $\varepsilon>0$. Then $F(n)\not\in
O\big{(}n^{(\log_{k}\rho)-\varepsilon}\big{)}$, $F(n)\in
O\big{(}n^{(\log_{k}\rho)+\varepsilon}\big{)}$ and $F_{0}^{s}\in
O\big{(}(k^{s})^{(\log_{k}\rho)+\varepsilon}\big{)}$ iff there is an
$S\in\mathbb{N}_{0}$ such that for all $s\geq S$, the inequalities

$(\log_{k}\rho)-\varepsilon<\log_{k}\rho_{s}\leq(\log_{k}\rho)+\varepsilon$

hold. But this is equivalent to
$\log_{k}\rho=\lim_{s\to\infty}\log_{k}\rho_{s}$ and therefore equivalent to

$\rho=\lim_{s\to\infty}\rho_{s},$

which completes the proof of the equivalence.

If $f(n)$ is as in the proposition, then by Proposition 3.5 we have the
equivalence of

$f(n)\in O(n^{\sigma}).$

to

$\lVert F(n)\rVert\in O(n^{\sigma})=O(k^{s\sigma})$

as $n\to\infty$ for any fixed real algebraic $\sigma$, so it is allowed to
replace $F(n)$ by $f(n)$ in our statements. ∎

###### Proof of Corollary G.

As $F_{0}$ is the zero-matrix, $F_{0}^{s}\in
O\big{(}(k^{s})^{(\log_{k}\rho^{\prime})+\varepsilon}\big{)}$ as $s\to\infty$
holds trivially, and the result therefore follows from Proposition 3.6. ∎

###### Proof of Theorem F.

We reduce from the recursively unsolvable question whether the joint spectral
radius of a set of matrices over the rationals is bounded by $1$; see Blondel
and Tsitsiklis [4].

So let us assume the decision problem of Theorem F is recursively solvable.
Let $F_{1}$, …, $F_{k-1}$ be square matrices over the rationals, let $F_{0}$
be the zero-matrix, and set $F(n)=F_{n_{0}}\cdots F_{n_{s-1}}$ for
$(n)_{k}=n_{s-1}\cdots n_{0}$. By the reduction of Lemma 3.4, where property
$P$ is whether a sequence does not have at least polynomial growth, we can
decide whether $F(n)$ does not have at least polynomial growth. This is
equivalent to deciding whether for all $\sigma>0$ we have $F(n)\in
O(n^{\sigma})$. By Corollary G with $\rho^{\prime}=1$, this is equivalent to
deciding whether the joint spectral radius $\rho$ of $F_{1}$, …, $F_{k-1}$
being at most $1$, a contradiction. ∎

##  4\. Images and preimages

By [1, Theorem 5.2], it is undecidable whether a given $k$-regular sequence
$(f(n))_{n\geq 0}$ has a zero term, i.e., whether there exists an
$n\in\mathbb{N}_{0}$ with $f(n)=0$.

###  4.1. Preimages

In this section we use closure properties of regular languages without further
comment. See, for example, [12].

######  Theorem H.

Let $k\geq 2$ be an integer and $\mathbb{S}\supseteq\mathbb{N}_{0}$ be a set.
The decision problem

> “Given a $k$-regular sequence $(f(n))_{n\geq 0}$ over $\mathbb{S}$ and a
> number $q\in\mathbb{S}$, is the language associated with $f^{-1}(q)$
> regular, i.e., can it be recognized by a deterministic finite automaton?”

is recursively unsolvable.

Above, the language associated with $f^{-1}(q)$ is
$\\{(n)_{k}\colon\mathopen{}f(n)=q\\}$, where $(n)_{k}$ is the standard
$k$-ary expansion of $n$.

A result from [1, Theorem 5.3] states that there exists a $k$-regular sequence
$(f(n))_{n\geq 0}$ such that neither $\\{(n)_{k}\colon\mathopen{}f(n)=0\\}$
nor $\\{(n)_{k}\colon\mathopen{}f(n)\neq 0\\}$ are context-free.

###### Proof of Theorem H.

We can assume that $q=0$ by subtracting $q$ from the $k$-regular sequence. In
order to prove the theorem, we reduce from Hilbert’s tenth problem.

For a given multivariate polynomial $p$ in $t$ variables over $\mathbb{Z}$, we
choose $r\in\mathbb{N}$ such that $K=k^{r}>t+2$, and we construct

(4.1) $f(n)=\bigl{(}\mathop{{#1}{}}(\lvert z\rvert_{1},\lvert
z\rvert_{2},\dots,\lvert z\rvert_{t})\bigr{)}^{2}+\bigl{(}\lvert
z\rvert_{t+1}-\lvert z\rvert_{t+2}\bigr{)}^{2},$

where $z=(n)_{K}$. The sequence $(f(n))_{n\geq 0}$ is $K$-regular (see Remark
1.1) and therefore $k$-regular as well by [1, Theorem 2.9].

The following claim shows that the above indeed provides a reduction.

######  Claim.

The set $f^{-1}(0)$ is not recognized by a deterministic finite automaton iff
there exist nonnegative integers $x_{1}$, $x_{2}$, …, $x_{t}$ such that
$p(x_{1},\dots,x_{t})=0$.

If $p(x_{1},\dots,x_{t})=0$ has no solution in $\mathbb{N}_{0}^{t}$, then
$f(n)\neq 0$ by its construction. Thus $f^{-1}(0)$ is the empty set and is
accepted by a deterministic finite automaton.

Otherwise, suppose we have nonnegative integers $x_{1}$, …, $x_{t}$ with
$p(x_{1},\dots,x_{t})=0$, and suppose the language
$L=\\{(n)_{K}\colon\mathopen{}f(n)=0\\}$ is accepted by a deterministic finite
automaton, i.e., $L$ is regular. We note that each $z=(n)_{K}\in L$ satisfies
$\lvert z\rvert_{t+1}=\lvert z\rvert_{t+2}$ as $f(n)=0$ and this is equivalent
to both squares in its definition (4.1) being zero. Moreover, for each
$s\in\mathbb{N}_{0}$, there is a $z\in L$ with $s=\lvert z\rvert_{t+1}=\lvert
z\rvert_{t+2}$.

As $L$ is regular, so is

$L_{1}=L\cap 1^{x_{1}}2^{x_{2}}\cdots t^{x_{t}}(t+1)^{+}(t+2)^{+}$

where $d^{+}=\\{d,dd,ddd,\ldots\\}$ for a letter $d$. Furthermore, the left
quotient

$L_{2}=\\{1^{x_{1}}2^{x_{2}}\cdots t^{x_{t}}(t+1)\\}^{-1}L_{1}$

is regular. However,

$L_{2}=\\{(t+1)^{m-1}(t+2)^{m}\colon\mathopen{}m\geq 1\\},$

which is not regular; see [7, Examples 2.8 and 5.2]. This contradiction proves
the desired result.

Now, if we can decide whether the language associated with $f^{-1}(q)$ is
regular, then we can decide whether a solution of $p$ exists, and therefore
decide Hilbert’s tenth problem. This completes the proof of Theorem H. ∎

###  4.2. Images

######  Theorem I.

Let $k\geq 2$ be an integer and $\mathbb{S}\supseteq\mathbb{N}_{0}$ be a set.
The decision problem

> “Given a $k$-regular sequence $(f(n))_{n\geq 0}$ over $\mathbb{S}$, is
> $\\{f(n)\colon\mathopen{}n\in\mathbb{N}_{0}\\}=\mathbb{N}_{0}$?”

is recursively unsolvable.

###### Proof.

In order to prove the theorem, we again reduce from Hilbert’s tenth problem.
For a given multivariate polynomial $p$ in $t$ variables over $\mathbb{Z}$, we
choose $r\in\mathbb{N}$ such that $K=k^{r}>t$, and we construct

$f(n)=\begin{cases}n/2+1,&\text{if $n$ is even;}\\\
\bigl{(}\mathop{{#1}{}}(\lvert z\rvert_{1},\lvert z\rvert_{2},\dots,\lvert
z\rvert_{t})\bigr{)}^{2},&\text{if $n$ is odd},\end{cases}$

where $z=\bigl{(}(n-1)/2\bigr{)}_{K}$. The sequence $(f(n))_{n\geq 0}$ is
$K$-regular (Remark 1.1 and [1, Theorem 2.7]) and therefore $k$-regular as
well by [1, Theorem 2.9].

Once we have shown the following claim, we have a reduction to Hilbert’s tenth
problem and therefore the proof of Theorem I is completed.

######  Claim.

The set $\\{f(n)\colon\mathopen{}n\in\mathbb{N}_{0}\\}$ equals
$\mathbb{N}_{0}$ iff there exist nonnegative integers $x_{1}$, $x_{2}$, …,
$x_{t}$ such that $p(x_{1},\dots,x_{t})=0$.

If $p(x_{1},\dots,x_{t})=0$ has no solution in $\mathbb{N}_{0}^{t}$, then
$f(n)\neq 0$ by its construction. Thus $0$ is not in the set
$\\{f(n)\colon\mathopen{}n\in\mathbb{N}_{0}\\}$, so this set cannot be equal
to $\mathbb{N}_{0}$.

Otherwise, suppose we have nonnegative integers $x_{1},\ldots,x_{t}$ with
$p(x_{1},\dots,x_{t})=0$, then there exists an $n\in\mathbb{N}_{0}$ with
$f(n)=0$. As $\\{f(n)\colon\mathopen{}\text{$n\in\mathbb{N}_{0}$ is even}\\}$
already contains all the positive integers, all nonnegative integers appear as
a value $f(n)$ somewhere. ∎

######  Theorem J.

Let $k\geq 2$ be an integer and $\mathbb{S}\supseteq\mathbb{Z}$ be a set. The
decision problem

> “Given a $k$-regular sequence $(f(n))_{n\geq 0}$ over $\mathbb{S}$, is
> $\\{f(n)\colon\mathopen{}n\in\mathbb{N}_{0}\\}=\mathbb{Z}$?”

is recursively unsolvable.

###### Proof.

The proof runs along the same lines as the proof of Theorem I, but for a given
multivariate polynomial $p$ in $t$ variables over $\mathbb{Z}$, we choose
$r\in\mathbb{N}$ such that $K=k^{r}>t$, and we construct

$f(n)=\begin{cases}n/3+1,&\text{if $n\equiv 0\pmod{3}$;}\\\
-(n-1)/3-1,&\text{if $n\equiv 1\pmod{3}$;}\\\ \mathop{{#1}{}}(\lvert
z\rvert_{1},\lvert z\rvert_{2},\dots,\lvert z\rvert_{t}),&\text{if $n\equiv
2\pmod{3}$,}\end{cases}$

where $z=\bigl{(}(n-2)/3\bigr{)}_{K}$. Then the set
$\\{f(n)\colon\mathopen{}n\in\mathbb{N}_{0}\\}$ equals $\mathbb{Z}$ iff there
exist nonnegative integers $x_{1}$, $x_{2}$, …, $x_{t}$ such that
$p(x_{1},\dots,x_{t})=0$. ∎

We can extend the above to the question whether two $k$-regular sequences have
the same image.

######  Corollary K.

Let $k\geq 2$ be an integer and $\mathbb{S}\supseteq\mathbb{N}_{0}$ be a set.
The decision problem

> “Given two $k$-regular sequences $(f(n))_{n\geq 0}$ and $(g(n))_{n\geq 0}$
> over $\mathbb{S}$, do their images coincide, i.e., is
> $\\{f(n)\colon\mathopen{}n\in\mathbb{N}_{0}\\}=\\{g(n)\colon\mathopen{}n\in\mathbb{N}_{0}\\}$?”

is recursively unsolvable.

###### Proof.

We reduce from the decision problem of Theorem I, so let $(f(n))_{n\geq 0}$ be
a $k$-regular sequence over $\mathbb{S}$ and set $g(n)=n$. If we can decide
whether these two sequences have the same image, then we decide whether

$\\{f(n)\colon\mathopen{}n\in\mathbb{N}_{0}\\}=\\{g(n)\colon\mathopen{}n\in\mathbb{N}_{0}\\}=\mathbb{N}_{0},$

which contradicts Theorem I. ∎

######  Theorem L.

Let $k\geq 2$ be an integer and $\mathbb{S}\supseteq\mathbb{N}_{0}$ be a set.
The decision problem

> “Given a $k$-regular sequence $(f(n))_{n\geq 0}$ over $\mathbb{S}$, does
> $f(n)$ take the same value twice?”

is recursively unsolvable.

###### Proof.

In order to prove the theorem, we reduce from Hilbert’s tenth problem. For a
given multivariate polynomial $p$ in $t$ variables over $\mathbb{Z}$, we
choose $r\in\mathbb{N}$ such that $K=k^{r}>t$, and we construct

$g(m)=\bigl{(}\mathop{{#1}{}}(\lvert z\rvert_{1},\lvert
z\rvert_{2},\dots,\lvert z\rvert_{t})\bigr{)}^{2}$

where $z=(m)_{K}$ and

$f(n)=\sum_{0\leq m<n}g(m)$

The sequence $(f(n))_{n\geq 0}$ is $K$-regular (Remark 1.1 and [1, Theorem
3.1]) and therefore $k$-regular as well by [1, Theorem 2.9].

Once we have shown the following claim, we have a reduction to Hilbert’s tenth
problem and therefore the proof of Theorem L is completed.

######  Claim.

The sequence $(f(n))_{n\geq 0}$ takes the same value twice iff there exist
nonnegative integers $x_{1}$, $x_{2}$, …, $x_{t}$ such that
$p(x_{1},\dots,x_{t})=0$.

If $p(x_{1},\dots,x_{t})=0$ has no solution in $\mathbb{N}_{0}^{t}$, then
$g(m)$ is strictly positive, and therefore $f(n)$ strictly increasing. So no
value is taken twice.

Otherwise, suppose we have nonnegative integers $x_{1}$, …, $x_{t}$ with
$p(x_{1},\dots,x_{t})=0$, then there exists an $n\in\mathbb{N}_{0}$ with
$g(n)=0$, and so $f(n)=f(n+1)$. ∎

##  5\. Squares and other $\alpha$-powers

Given a sequence $(f(n))_{n\geq 0}$ and an integer $\alpha\geq 2$, an
$\alpha$-power is a nonempty contiguous subsequence $(f(j))_{i\leq j<i+\alpha
m}$ of length $\alpha$, for some $i$ and $m$, such that $f(i+t)=f(i+sm+t)$ for
all $0\leq s<\alpha$ and $0\leq t<m$. We call a $2$-power a square. For
example, the fractional part of the decimal representation of $e$ contains the
square $18281828$.

A palindrome is a nonempty contiguous subsequence that reads the same forwards
and backwards. A palindrome is nontrivial if it is of length $\geq 2$.

For automatic sequences, the presence of squares, higher powers, and
nontrivial palindromes is decidable (see, e.g., [8]). We now show that, in
contrast, the existence of these patterns is undecidable for $k$-regular
sequences.

######  Theorem M.

Let $\alpha\geq 2$ be an integer. The decision problem

> “Given a $k$-regular sequence $(f(n))_{n\geq 0}$, does $(f(n))_{n\geq 0}$
> contain an $\alpha$-power?”

is recursively unsolvable.

###### Proof.

We reduce from the problem of deciding whether a $k$-regular sequence has a
$0$ term.

Given a $k$-regular sequence $(f(n))_{n\geq 0}$ for which we want to decide
whether $f(n)=0$ for some $n$, we can (effectively) transform it to the
$k$-regular sequence $(g(n))_{n\geq 0}$ defined recursively by $g(0)=1$ and

$g(n)=g(n-1)+f(n-(\alpha-1))^{2}\cdots f(n-2)^{2}\,f(n-1)^{2},\quad\text{for
$n\geq 1$}.$

(Note that we use the convention $f(-i)=1$ for $i\geq 1$.) For squares, this
simplifies to the explicit formula $g(n)=1+f(0)^{2}+\cdots+f(n-1)^{2}$. Then
$(g(n))_{n\geq 0}$ is (not necessarily strictly) increasing, so it contains an
$\alpha$-power iff there exists $n\geq 0$ such that
$g(n)=g(n+1)=\cdots=g(n+\alpha-1)$. But this occurs iff $f(n)$ = 0. ∎

Using the same technique, we can prove following theorem for palindromes.

######  Theorem N.

The decision problem

> “Given a $k$-regular sequence $(f(n))_{n\geq 0}$, does $(f(n))_{n\geq 0}$
> contain a nontrivial palindrome?”

is recursively unsolvable.

###### Proof.

The same proof given for squares above works unchanged. ∎

## References

  * [1] Jean-Paul Allouche and Jeffrey Shallit, _The ring of $k$-regular sequences_, Theoret. Comput. Sci. 98 (1992), no. 2, 163–197.
  * [2] Jean-Paul Allouche and Jeffrey Shallit, _Automatic sequences: Theory, applications, generalizations_ , Cambridge University Press, 2003.
  * [3] Jean Berstel and Christophe Reutenauer, _Noncommutative rational series with applications_ , Encyclopedia of Mathematics and its Applications, vol. 137, Cambridge University Press, Cambridge, 2011.
  * [4] Vincent D. Blondel and John N. Tsitsiklis, _The boundedness of all products of a pair of matrices is undecidable_, Systems Control Lett. 41 (2000), no. 2, 135–140.
  * [5] Alan Cobham, _Uniform tag sequences_ , Math. Systems Theory 6 (1972), 164–192.
  * [6] Martin Davis, Hilary Putnam, and Julia Robinson, _The decision problem for exponential diophantine equations_, Ann. of Math. (2) 74 (1961), 425–436.
  * [7] Samuel Eilenberg, _Automata, languages, and machines. Vol. A_ , Academic Press [A subsidiary of Harcourt Brace Jovanovich, Publishers], New York, 1974, Pure and Applied Mathematics, Vol. 58.
  * [8] Émilie Charlier, Narad Rampersad, and Jeffrey Shallit, _Enumeration and decidable properties of automatic sequences_ , Internat. J. Found. Comp. Sci. 23 (2012), 1035–1066.
  * [9] Albrecht Fröhlich and John C. Shepherdson, _Effective procedures in field theory_, Philos. Trans. Roy. Soc. London Ser. A 248 (1956), 407–432.
  * [10] Edward R. Griffor (ed.), _Handbook of computability theory_ , Studies in Logic and the Foundations of Mathematics, vol. 140, North-Holland Publishing Co., Amsterdam, 1999.
  * [11] Juha Honkala, _Quasi-universal $k$-regular sequences_, Theoret. Comput. Sci. 891 (2021), 84–89.
  * [12] John E. Hopcroft and Jeffrey D. Ullman, _Introduction to automata theory, languages, and computation_ , Addison-Wesley, 1979.
  * [13] Gérard Jacob, _Décidabilité de la finitude des demi-groupes de matrices_ , Theoretical computer science (Third GI Conf., Darmstadt, 1977), Lecture Notes in Computer Science, vol. 48, Springer, Berlin, 1977, pp. 259–269.
  * [14] by same author, _Un algorithme calculant le cardinal, fini ou infini, des demi-groupes de matrices_, Theoret. Comput. Sci. 5 (1977/78), no. 2, 183–204.
  * [15] by same author, _La finitude des représentations linéaires des semi-groupes est décidable_, J. Algebra 52 (1978), no. 2, 437–459.
  * [16] Arnaldo Mandel and Imre Simon, _On finite semigroups of matrices_, Theoret. Comput. Sci. 5 (1977/78), no. 2, 101–111.
  * [17] Yuri V. Matijasevič, _The Diophantineness of enumerable sets_ , Dokl. Akad. Nauk SSSR 191 (1970), 279–282.


# Detection of the extended $\gamma$-ray emission around supernova remnant DA
530 with Fermi-LAT

Yuliang Xin Xiaolei Guo School of Physical Science and Technology, Southwest
Jiaotong University, Chengdu 610031, China
<EMAIL_ADDRESS>
###### Abstract

We report the extended GeV $\gamma$-ray emission around the high Galactic
latitude supernova remnant (SNR) DA 530 with the PASS 8 data recorded by the
Fermi Large Area Telescope (Fermi-LAT). The $\gamma$-ray spectrum in the
energy range of 100 MeV - 1 TeV follows a power law model with an index of
2.23. The much more extended $\gamma$-ray emission than the radio shell of DA
530 and the spatial coincidence with the molecular cloud suggest that the
$\gamma$-ray emission could be originated from the hadronic process, where the
high energy protons are accelerated in and escaped from the shock of DA 530.
With a steady-state injection model of protons, the $\gamma$-ray spectrum can
be well fitted with the typical Galactic value for diffusion coefficient and
the low energy content of the total escaped protons.

gamma rays: general - gamma rays: ISM - ISM: individual objects (DA 530) -
radiation mechanisms: non-thermal

##  1 Introduction

Supernova remnants (SNRs) are widely believed to be the dominant accelerators
of Galactic cosmic rays (CRs). And CRs can be accelerated by the high speed
shock of SNR with the mechanism of diffusive shock acceleration (DSA;
Blandford & Eichler, 1987). The highest-energy CRs are expected to escape from
the shock due to the absence of the self-generated magnetic turbulence
(Ptuskin & Zirakashvili, 2005; Fujita et al., 2011). If there are some
molecular clouds (MCs) in the vicinity of SNRs, these MCs could be illuminated
by the escaped protons and produce the intense $\gamma$-ray emission with the
hadronic process, i.e., the $\gamma$-ray emission are believed to be from the
decay of neutral pions produced in inelastic collisions between escaped
protons and the dense gas in MCs. And several sources have been detected in
this perspective, e.g. W28 (Aharonian et al., 2008; Cui et al., 2018), W44
(Uchiyama et al., 2012), SNR G15.4+0.1 (Li et al., 2023), SNR G45.7-0.4 (Zhang
et al., 2021), etc. Especially in the case of W28, the GeV and TeV
$\gamma$-ray emission are detected in three distinct regions, which are
spatially coinciding with MCs offset from the shock of SNR (Aharonian et al.,
2008; Cui et al., 2018). Searching for the $\gamma$-ray emission from such
associations could be helpful to explore the propagation properties of high
energy particles escaped from SNRs and afterwards travelled through the
interstellar medium (ISM), which further provides a key ingredient to
establishing the connection between SNRs and the origin of Galactic CRs.

DA 530, also known as G93.3+6.9, is a high Galactic latitude SNR (Roger &
Costain, 1976). And based on the radio observations, DA 530 is classified to
be a shell-type SNR with a bilateral morphology (Landecker et al., 1999). The
radio radiation of DA 530 has extremely high polarization percentage, reaching
more than 50%, which could be interpreted by the well-ordered magnetic field
across the whole remnant (Haslam et al., 1980; Lalitha et al., 1984). The
X-ray emission of DA 530 first detected by ROSAT is extremely faint with a
centrally brightened morphology (Landecker et al., 1999). By re-analysing the
Chandra data, Jiang et al. (2007) found a small-scale hard X-ray feature near
the centre of the remnant, which is argued to be a pulsar wind nebula (PWN)
associated with SNR. And the age of DA 530 was suggested to be $\sim$ 5000 yrs
based on the canonical blast wave model from Sedov (1959). However, the
subsequent XMM-Newton observations detected a large extended source (XMM
J205314.4+551528) in the radio bright southeast (SE) rim of DA 530 (Bocchino
et al., 2008), and the authors explained it to be the PWN associated with DA
530. The recent Suzaku data analysis of the SE rim of DA 530 confirmed the
results of XMM-Newton, which suggested that the PWN scenario can not be ruled
out (Deniz et al., 2022). Nonetheless, the searching for the pulsars or
compact central sources associated with DA 530 in the radio and X-ray bands
show the null results (Lorimer et al., 1998; Landecker et al., 1999; Kaplan et
al., 2004; Straal & van Leeuwen, 2019). The distance of DA 530 is not very
clear so far, which was first derived to be 6.9 $\pm$ 2.2 kpc using the
empirical surface brightness-diameter ($\Sigma$ \- D) relation of SNRs (Roger
& Costain, 1976). Then an updated distance of 2 - 5 kpc was given by Haslam et
al. (1980). Based on the neutral hydrogen (HI) observations with Dominion
Radio Astrophysical Observatory Synthesis Telescope (DRAO-ST), Landecker et
al. (1999) reported that DA 530 lies within a shell of HI, possibly created by
an earlier stellar wind of the progenitor. And the distance of it is estimated
to be 1.0 - 3.5 kpc. Subsequently, Foster & Routledge (2003) derived a
distance of 2.2 $\pm$ 0.5 kpc using the updated method for the absorption
column density. Using the data from DRAO-ST and National Radio Astronomy
Observatory Very Large Array (NRAO-VLA), Booth et al. (2022) observed the
absorption by intervening HI of the polarized emission from DA 530, and
concluded that the distance of DA 530 is 4.4${}^{\rm+0.4}_{\rm-0.2}$ kpc.
Using the Seoul Radio Astronomy Observatory (SRAO) 6-m telescope CO
observations, Jeong et al. (2012) detected CO emission at -6 to +5 km s-1 at
the northeast boundary of DA 530, which shows that there is a large diffuse
molecular cloud.

In this work, we report the detection of the extended $\gamma$-ray emission
around SNR DA 530, with the PASS 8 data recorded by Fermi-LAT. The data
analysis method and results are shown in Section 2, including the spatial and
spectral analyses. The observations of molecular cloud around DA 530 is
presented in Section 3. And the discussion of the potential origin of the
$\gamma$-ray emission is shown in Section 4, followed by the summary in
Section 5.

##  2 Fermi-LAT Data Analysis

###  2.1 Data Reduction

Fermi-LAT is a pair-conversion $\gamma$-ray telescope that is sensitive to
photon energies greater than 20 MeV. The LAT has continuously monitored the
sky since 2008 and scans the entire sky every 3 hr (Atwood et al., 2009). Note
that the latest released Pass 8 data set has significant improvements in
comparison with the former ones, including an enhanced effective area,
especially in the low energy range and better point-spread function
(PSF)111https://www.slac.stanford.edu/exp/glast/groups/canda/lat$\\_$Performance.htm.
And in the following analysis, we select the latest Pass 8 version of Fermi-
LAT data recorded from August 4, 2008 (Mission Elapsed Time 239557418) to
August 4, 2022 (Mission Elapsed Time 681264005) with “Source” event class
(evclass = 128 & evtype = 3) to analyse the $\gamma$-ray emission around DA
530. The region of interest (ROI) is a $20^{\circ}\times 20^{\circ}$ square
region centered at the position of DA 530 (R.A. = $313^{\circ}\\!.14$, decl. =
$55^{\circ}\\!.36$; Roger & Costain, 1976). And in order to reduce the
contamination from Earth Limb, the events with zenith angle larger than
$90^{\circ}$ are excluded. We adopt the events with energy range of 100 MeV -
1 TeV for the spectral analysis. While for the spatial analysis, the events in
the energy range of 1 GeV - 1 TeV is selected considering the impact of PSF of
Fermi-LAT. The data are analyzed using the standard Fermi ScienceTools
222http://fermi.gsfc.nasa.gov/ssc/data/analysis/software/ with the
instrumental response function (IRF) of “P8R3_SOURCE_V3”. The binned
likelihood analysis method with gtlike is used to fit the data. To model the
Galactic and isotropic diffuse background emissions, gll_iem_v07.fits and
iso_P8R3_SOURCE_V3_v1.txt
333http://fermi.gsfc.nasa.gov/ssc/data/access/lat/BackgroundModels.html are
adopted. All sources in the incremental version of the fourth Fermi-LAT source
catalog (4FGL-DR3; Abdollahi et al., 2020, 2022) within a radius of
$20^{\circ}$ from the ROI center and the two components of the diffuse
background, are included in the source model, which is generated by the user-
contributed software
make4FGLxml.py444http://fermi.gsfc.nasa.gov/ssc/data/analysis/user/. During
the likelihood analysis, the normalizations and the spectral parameters of all
sources within $7^{\circ}$ to the center of ROI, together with the
normalizations of the two components of the diffuse background, are set to be
free.

###  2.2 Spatial Analysis

In the region of DA 530, a $\gamma$-ray point source (4FGL J2051.1+5539) is
listed in the 4FGL-DR3 catalog, which has no identified counterpart (Abdollahi
et al., 2022). First, we create a $3^{\circ}\\!.0$ $\times$ $3^{\circ}\\!.0$
Test Statistic (TS) map by subtracting the emission from the sources and
backgrounds (except 4FGL J2051.1+5539) in the best-fit model with gttsmap,
which is shown in the left panel of Figure 1. At the center of ROI, a
significant $\gamma$-ray excess (labelled as SrcT) is found around SNR DA 530.
For the spatial template of SrcT, we first treat it to be a point source, and
the best-fit position is given to be R.A. = $312^{\circ}\\!.983\pm
0^{\circ}\\!.029$, decl. = $55^{\circ}\\!.666\pm 0^{\circ}\\!.029$ by adopting
Fermipy, a PYTHON package that automates analyses with Fermi ScienceTools
(Wood et al., 2017). The TS value of SrcT as a point source is 33.46 with the
new coordinate. Then we carried out a spatial extension test for the
$\gamma$-ray emission of SrcT using an uniform disk and a two-dimensional (2D)
Gaussian template with Fermipy. And the best-fit central positions and
extensions of the spatial templates are listed in Table 1, together with the
corresponding TS value of SrcT and the maximum likelihood values. We compared
the overall maximum likelihood of the extended template ($\rm\mathcal{L}_{\rm
ext}$) with that of the point source model ($\rm\mathcal{L}_{\rm pt}$), and
the significance of the extended model is defined to be $\rm{TS_{\rm
ext}}=2(\ln\mathcal{L}_{\rm ext}-\ln\mathcal{L}_{\rm pt})$. Lande et al.
(2012) suggests that a source can be assessed to be significantly extended if
$\rm{TS_{\rm ext}}>16$. The different maximum likelihood values of the
different templates show that an uniform disk can best fit the $\gamma$-ray
emission from SrcT. And the central position and the 68% containment radius of
the uniform disk are fitted to be R.A. = $313^{\circ}\\!.611$, decl. =
$55^{\circ}\\!.344$ and R68 = $0^{\circ}\\!.527$, respectively. The value of
$\rm{TS_{\rm ext}}$ between the uniform disk model and point source model is
calculated to be 26.2, corresponding to $\sim$5.1$\sigma$ extension with one
additional degree of freedom (dof). With the uniform disk template, the TS
value of SrcT is fitted to be 58.75 in the energy range of 1 GeV - 1 TeV,
corresponding to a significance level of $\sim$6.7 $\sigma$ with five degrees
of freedom.

Table 1: Spatial Analysis for SrcT in the energy range of 1 GeV - 1 TeV Spatial Model | R.A., decl. | TS Value | Degrees of Freedom | -log(Likelihood)
---|---|---|---|---
Point Source |  $312^{\circ}\\!.983\pm 0^{\circ}\\!.029$, $55^{\circ}\\!.666\pm 0^{\circ}\\!.029$ | 33.46 | 4 | -182258.99
Uniform disk |  $313^{\circ}\\!.611\pm 0^{\circ}\\!.043$, $55^{\circ}\\!.344\pm 0^{\circ}\\!.049$, R68 = $0^{\circ}\\!.527^{+0^{\circ}\\!.028}_{-0^{\circ}\\!.029}$ | 58.75 | 5 | -182272.08
2D-Gaussian |  $313^{\circ}\\!.499\pm 0^{\circ}\\!.100$, $55^{\circ}\\!.415\pm 0^{\circ}\\!.102$, R68 = $0^{\circ}\\!.567^{+0^{\circ}\\!.101}_{-0^{\circ}\\!.080}$ | 53.93 | 5 | -182269.65

Figure 1: Left: $3^{\circ}\\!.0$ $\times$ $3^{\circ}\\!.0$ TSmap in the energy
range of 1 GeV - 1 TeV. The red cross shows the position of 4FGL J2051.1+5539
in 4FGL-DR3 catalog. And the green solid circle marks the best-fit 68%
containment radius of the uniform disk for the spatial template of SrcT. The
radio image of SNR DA 530 at 1.4 GHz is shown as the cyan contours. Right: SED
of SrcT in the energy range of 100 MeV - 1 TeV with the gray histogram shown
as the TS value for each energy bin. The red error bars show the statistical
errors and the sums of the statistical and systematic errors calculated by
$\sigma$ = $\sqrt{\sigma_{\rm stat}^{2}+\sigma_{\rm syst}^{2}}$ are marked by
the blue error bars. The arrows indicate the 95% upper limits for the energy
bin with TS value of SrcT smaller than 5.0. The black solid and dashed lines
show the global best-fit power law spectrum and its 1$\sigma$ statistic error
in the energy range of 100 MeV - 1 TeV.

###  2.3 Spectral Analysis

To investigate the $\gamma$-ray spectrum of SrcT, the global likelihood
analysis is performed in the energy range from 100 MeV to 1 TeV with the
spatial template of an uniform disk. The spectrum of SrcT can be well
described by a power law (PL) model. The spectral index and the integral
photon flux in the energy range of 100 MeV - 1 TeV are fitted to be 2.23 $\pm$
0.09 and $(1.25\pm 0.23)\times 10^{-8}$ photon cm-2 s-1, respectively. And we
also adopted a log-parabola model (LPb) to test the spectral curvature of
SrcT. The variation of the maximum likelihood values between PL and LPb models
is only $\rm{TS_{\rm curve}}=2(\ln\mathcal{L}_{\rm LPb}-\ln\mathcal{L}_{\rm
PL})=2.1$, which suggests no significant curvature for the $\gamma$-ray
spectrum of SrcT. With the PL model for SrcT, we divide the data into ten
logarithmically equal energy bins from 100 MeV to 1 TeV and repeat the
likelihood fitting to give the spectral energy distribution (SED) of SrcT. For
the likelihood analysis, only the spectral normalizations of sources within
$7^{\circ}$ from SrcT are left free, together with the normalizations of the
two components of the diffuse background. And the spectral indices of these
sources are fixed to be the best-fit values in the global likelihood analysis.
For the energy bin with a TS value of SrcT smaller than 5.0, an upper limit
with a 95% confidence level is calculated. For the SED, we also estimated the
systematic errors due to the Galactic diffuse emission by changing the best-
fit normalization of the Galactic diffuse model artificially by $\pm$6% (Abdo
et al., 2010). And the sums of the statistical and systematic errors are
calculated by $\sigma$ = $\sqrt{\sigma_{\rm stat}^{2}+\sigma_{\rm syst}^{2}}$
for each energy bin. The SED of SrcT is shown in the right panel of Figure 1,
which is also consistent with the global fitting of the power law model.

##  3 CO Observation

With the CO observations from SRAO, Jeong et al. (2012) claimed a large
diffuse molecular cloud surrounding the northeast boundary of DA 530. However,
the size of the $\gamma$-ray emission detected here is much larger than that
of their CO observation. Therefore, we try to search for the components of
molecular cloud in a large extend using the data from the CfA 1.2m millimeter-
wave telescope to understand the origin of the $\gamma$-ray excess (Dame et
al., 2001). And we found that in the $\gamma$-ray emission region, the
velocity distribution of the CO content shows a clear excess in the velocity
range of -6 $\sim$ +5 km s-1 as shown in Figure 2, which is consistent with
the range in Jeong et al. (2012). Adopting the standard Galactic rotation
model (Reid et al., 2016, 2019), the velocity interval corresponds to a
kinetic distance of $\sim$1.6 kpc. And considering the systematic
uncertainties due to the rotation curve, the derived distance is close to the
value of 2.2 $\pm$ 0.5 kpc derived by Foster & Routledge (2003) using the
updated method for the absorption column density, which is also much lower
that the result of 4.4 kpc in Booth et al. (2022).

We have then estimated the mass content of the molecular material in the
$\gamma$-ray emission region with

$M=\mu
m_{\mathrm{H}}d^{2}\Omega_{\mathrm{px}}X_{\mathrm{CO}}{\sum_{\mathrm{px}}}W_{\mathrm{CO}}$
(1)

Here the mean molecular weight $\mu$ is adopted to be 2.8 assuming a relative
helium abundance of 25%. $m_{\mathrm{H}}$ is the mass of the Hydrogen nucleon,
and the distance is adopted to be d = 2.2 kpc (Foster & Routledge, 2003).
$\Omega_{\mathrm{px}}$ is the solid angle subtended for each pixel in the map
shown as in Figure 2. And the value of the conversion factor of
$X_{\mathrm{CO}}=2\times 10^{20}$ cm-2 (K km s${}^{-1})^{-1}$ is adopted here
(Bolatto et al., 2013). $\sum_{\mathrm{px}}$WCO is calculated by summing the
map content of each pixel in the desired sky region and velocity range. For
the region of $0.527\arcdeg$ sky integration radius, the total mass of
molecular clouds in the region within the 68% containment radius of the
extended $\gamma$-ray source is estimated to be $\sim$ $2.6\times
10^{5}d_{2.2}^{2}M_{\odot}$. And the corresponding average gas number density
is about $\rm{n_{\rm gas}}$ = 300 cm-3 by assuming a spherical geometry of the
gas distribution.

Figure 2: Integrated CO emission intensity (K km s-1) around DA 530 in the
velocity range of -6 $\sim$ +5 km s-1. The 68% containment radius of the
extended $\gamma$-ray emission of SrcT is marked by the yellow solid circle.
The green and cyan contours show the $\gamma$-ray morphology of SrcT and the
radio image of SNR DA 530 at 1.4 GHz, as shown in the left panel of Figure 1.

##  4 Discussion

The Fermi-LAT data analysis above shows an extended $\gamma$-ray source SrcT,
which locates around the high Galactic latitude SNR DA 530. And the GeV
$\gamma$-ray spectrum of SrcT can be described by a power law model with an
index of 2.23 $\pm$ 0.09. To determine the origin of the $\gamma$-ray
emission, we also searched for the CO observation, and SrcT is spatially
consistent with the molecular clouds, which is located outside the radio shell
of DA 530. The spatial coincidence between the $\gamma$-ray emission and the
molecular clouds suggests that the $\gamma$-ray emission from SrcT could be
produced by the hadronic $\pi^{0}$ decay originating from the interaction
between the molecular gas and high energy protons, which are accelerated in
and escaped from the shock of DA 530. Such scenario is similar to the origin
of the $\gamma$-ray emission around SNR G15.4+1.0 (Li et al., 2023) and the
SNR associated with PSR J0837-2454(Zhang & Xin, 2023).

To explain the $\gamma$-ray emission from SrcT, we assume the steady-state
injection of protons into an uniform emission region. And the injection time
is adopted to be the age of DA 530 with T = 5000 yrs from Jiang et al. (2007).
The spectrum of injected protons is adopted to be a power law with an
exponential cutoff:

$\displaystyle Q(E)={Q_{0}}E^{-\Gamma}\exp\left(-\frac{E}{E_{\rm
p,cut}}\right).$ (2)

Here the spectral index is suggested to be $\Gamma=2.0$, which is consistent
with the radio spectral index of DA 530 of $\alpha$ = 0.45 $\pm$ 0.04 by
assuming the same spectral index for electrons and protons accelerated by the
shock of SNR. The cutoff energy of protons can not be well constrained and is
first adopted to be the energy of the cosmic ray knee with $E_{\rm p,cut}$ = 3
PeV. And the total energy of the injected protons is assumed to be $\rm W_{\rm
p,inj}$ = $\rm\eta E_{\rm SN}$, where $\eta$ is the fraction of the kinetic
energy of DA 530 converted into the escaped proton energy, and the the kinetic
energy of DA 530, $\rm E_{\rm SN}$, is adopted to be a typical value of 1051
erg.

By integrating the Eq. 16 in Thoudam & Hörandel (2012) on the variable radius
over a sphere with radius R, the escaped proton spectrum within the
$\gamma$-ray emission region can be derived as (Aharonian & Atoyan, 1996;
Thoudam & Hörandel, 2012):

$N_{p}(E,T)=\frac{Q(E)}{4\pi D(E)T}\int_{0}^{R}4\pi rdr\;{\rm
erfc}\left[\frac{r}{\sqrt{4D(E)T}}\right]$ (3)

Here, the diffusion coefficient of protons is set to be spatially constant and
energy-dependent with $D(E)=\chi D_{0}(E/E_{0})^{\delta}$, where
$D_{0}=3\times 10^{28}$ cm2 s-1 at $E_{0}=10$ GeV and $\chi$ = 1.0 corresponds
to the typical value of Galactic diffusion coefficient (Blasi, 2013). The
value of $\delta$ is adopted to be 1/3 or 1/2, which corresponds to the
Kolmogorov turbulence or Kraichnan turbulence for the diffusion coefficient
(Ptuskin et al., 2006). And with the distance of 2.2 kpc and 68% containment
radius of $0^{\circ}\\!.527$ for SrcT, the physical radius of the $\gamma$-ray
emission region is estimated to be R = 20.2 pc. For an injected spectrum given
by $Q(E)\propto E^{-\Gamma}$ and $D(E)\propto E^{\delta}$, the spectrum of
escaped protons, $N_{p}(E)$, approximately follows $N_{p}(E)\propto
E^{-(\Gamma+\delta)}$ in the high energy, where the diffusion radius of
protons defined as $\sqrt{4D(E)T}$ is larger than the size of the emission
region R. The typical values of $\eta$, $\chi$ and $\delta$ are selected to
calculate the different spectra of escaped protons in the $\gamma$-ray
emission region, and the corresponding $\gamma$-ray fluxes are calculated
using the $naima$ package by adopting the ambient gas density of $\rm{n_{\rm
gas}}$ = 300 cm-3 (Zabalza, 2015).

Figure 3 shows the resulting hadronic $\gamma$-ray spectra with the different
parameters. The spectra with the typical value of Galactic diffusion
coefficient, $\chi$ = 1.0, could explain the observational GeV data. And the
total energy of injected protons above 1 GeV is fitted to be about 2.0
$\times$ 1048 erg. which is lower than the estimated value by assuming the
total acceleration efficiency of 5% - 10% and the typical kinetic energy of
1051 erg for SNR (Blasi, 2013). Such result could be attributed to the
possible low kinetic energy of DA 530 with $\sim$1049 erg (Jiang et al.,
2007). And it also could mean that the bulk of accelerated particles are still
trapped inside the remnant, and the non-significant $\gamma$-ray emission
within the shell of DA 530 shown in the left panel of Figure 1 could be
explained by the low gas density inside SNR, which is also consistent with the
evolution environment of a stellar wind bubble for DA 530 (Landecker et al.,
1999). The total energies of escaped protons above 1 GeV in the $\gamma$-ray
emission region are calculated to be 1.1 and 1.3 $\times$ 1047 ($n_{\rm
gas}$/300 ${\rm cm^{-3}}$)-1 erg for $\delta$ = 1/2 and $\delta$ = 1/3,
respectively. To increase the total energy of injected protons, one needs to
increase the diffusion coefficient to make the $\gamma$-ray spectrum not be
changed significantly. By fixing the value of $\rm W_{\rm p,inj}$ = 1050 erg,
$\eta$ = 0.1, the diffusion coefficient needs to be about two orders of
magnitude higher than the typical Galactic value. And the estimated total
energy of escaped protons in the $\gamma$-ray emission region is about
1.2$\times$1047 ($n_{\rm gas}$/300 ${\rm cm^{-3}}$)-1 erg. In addition, we
also decreased the cutoff energy of protons to estimate the different escaping
models. And the allowed minimum value of cutoff energy is about 10 TeV, shown
as the red dotted line in Figure 3 with the total energy of escaped protons in
the $\gamma$-ray emission region of 1.1$\times$1047 ($n_{\rm gas}$/300 ${\rm
cm^{-3}}$)-1 erg, which expects a much lower flux in the TeV band.

Moreover, we also considered the the minimum energy of protons that can escape
from the SNR with a simple but reasonable approach of $E_{\rm esc}=E_{\rm
max}(t/t_{\rm sed})^{-\alpha}$ (Gabici et al., 2009; Ohira et al., 2011;
Thoudam & Hörandel, 2012), which assuming that the escape of the highest
energy particles starts at the onset of the Sedov phase of SNR. Here the
maximum energy of protons, $E_{\rm max}$, is assumed to be equal to $E_{\rm
p,cut}$ with 3 PeV, and $t_{\rm sed}$ = 500 yrs (Thoudam & Hörandel, 2012).
SNR DA 530 with the age of 5000 yrs is suggested to be in the Sedov phase
(Jiang et al., 2007), and the values of $E_{\rm esc}$ are calculated to be in
the range of 30 TeV - 2 PeV for the different values of $\alpha$ from 0.2 to
2.0 (Thoudam & Hörandel, 2012). And such energy range of escaping protons is
in accord with the needed one that explaining the GeV $\gamma$-ray spectrum.

Figure 3: Modeling of the $\gamma$-ray spectra with the hadronic escaping
models. The red solid, green and blue dashed lines indicate the scenarios with
the different values for $\chi$, $\delta$ and $\eta$ and the cutoff energy of
$E_{\rm p,cut}$ = 3 PeV. And the red dotted line shows the scenario with
$E_{\rm p,cut}$ = 10 TeV, as shown in the legend. The cyan dot–dashed line
shows the differential sensitivity of CTA-North (50 hr; Cherenkov Telescope
Array Consortium et al., 2019).

##  5 Summary

In this work, we analyzed the GeV $\gamma$-ray emission in the field of SNR DA
530, using 14 years of Fermi-LAT data, and found an extended $\gamma$-ray
source, SrcT, around DA 530, which can be described by an uniform disk
template. The GeV $\gamma$-ray spectrum of SrcT can be fitted by a power law
model with an index of 2.23 $\pm$ 0.09. The size of the $\gamma$-ray emission
region is much larger than that of the radio shell of DA 530. Based on the CO
observation from the CfA 1.2m millimeter-wave telescope, we found that the
molecular cloud component is spatially consistent with the $\gamma$-ray
emission region. Considering the much more extended $\gamma$-ray emission and
the spatial coincidence with the molecular cloud, the $\gamma$-ray emission of
SrcT is suggested to be from the hadronic $\pi^{0}$ decay due to the inelastic
collisions between the cloud and the high energy protons, which are
accelerated in and escaped from the shock of DA 530. With the assumption of
steady-state injection of protons, the $\gamma$-ray spectrum can be well
explained by the model with the typical Galactic value for the diffusion
coefficient ($\chi=1.0$). The total energy of the injected protons needs to be
much lower, which can be explained by the low kinetic energy of DA 530 or the
assumption that the bulk of accelerated particles are still trapped inside the
remnant. The hadronic $\gamma$-ray spectra with the different parameters
expect the different fluxes in the TeV band. And the potential detection by
the Cherenkov Telescope Array in the northern hemisphere (CTA-North; Cherenkov
Telescope Array Consortium et al., 2019) in the future could help to test the
different models and to constrain the escaping proton energy. Moreover, the
high-resolution observations for molecular cloud around DA 530 are necessary
to clearly identify the origin of the extended $\gamma$-ray emission.

We would like to thank the anonymous referee for very helpful comments, which
help to improve the paper. This work is supported by the National Natural
Science Foundation of China under the grants 12103040, 12147208 and U1931204,
and the Natural Science Foundation for Young Scholars of Sichuan Province,
China (No. 2022NSFSC1808).

## References

  * Abdo et al. (2010) Abdo, A. A., Ackermann, M., Ajello, M., et al. 2010, ApJ, 714, 927, doi: 10.1088/0004-637X/714/1/927
  * Abdollahi et al. (2020) Abdollahi, S., Acero, F., Ackermann, M., et al. 2020, ApJS, 247, 33, doi: 10.3847/1538-4365/ab6bcb
  * Abdollahi et al. (2022) Abdollahi, S., Acero, F., Baldini, L., et al. 2022, ApJS, 260, 53, doi: 10.3847/1538-4365/ac6751
  * Aharonian et al. (2008) Aharonian, F., Akhperjanian, A. G., Bazer-Bachi, A. R., et al. 2008, A&A, 481, 401, doi: 10.1051/0004-6361:20077765
  * Aharonian & Atoyan (1996) Aharonian, F. A., & Atoyan, A. M. 1996, A&A, 309, 917
  * Atwood et al. (2009) Atwood, W. B., Abdo, A. A., Ackermann, M., et al. 2009, ApJ, 697, 1071, doi: 10.1088/0004-637X/697/2/1071
  * Blandford & Eichler (1987) Blandford, R., & Eichler, D. 1987, Phys. Rep., 154, 1, doi: 10.1016/0370-1573(87)90134-7
  * Blasi (2013) Blasi, P. 2013, A&A Rev., 21, 70, doi: 10.1007/s00159-013-0070-7
  * Bocchino et al. (2008) Bocchino, F., Bandiera, R., & Gallant, Y. 2008, Advances in Space Research, 41, 407, doi: 10.1016/j.asr.2007.06.036
  * Bolatto et al. (2013) Bolatto, A. D., Wolfire, M., & Leroy, A. K. 2013, ARA&A, 51, 207, doi: 10.1146/annurev-astro-082812-140944
  * Booth et al. (2022) Booth, R. A., Kothes, R., Landecker, T., et al. 2022, ApJ, 941, 17, doi: 10.3847/1538-4357/ac9cd9
  * Cherenkov Telescope Array Consortium et al. (2019) Cherenkov Telescope Array Consortium, Acharya, B. S., Agudo, I., et al. 2019, Science with the Cherenkov Telescope Array, doi: 10.1142/10986
  * Cui et al. (2018) Cui, Y., Yeung, P. K. H., Tam, P. H. T., & Pühlhofer, G. 2018, ApJ, 860, 69, doi: 10.3847/1538-4357/aac37b
  * Dame et al. (2001) Dame, T. M., Hartmann, D., & Thaddeus, P. 2001, ApJ, 547, 792, doi: 10.1086/318388
  * Deniz et al. (2022) Deniz, C., Sezer, A., Bakis, H., & Raycheva, N. C. 2022, Advances in Space Research, 69, 2342, doi: 10.1016/j.asr.2021.12.024
  * Foster & Routledge (2003) Foster, T., & Routledge, D. 2003, ApJ, 598, 1005, doi: 10.1086/378947
  * Fujita et al. (2011) Fujita, Y., Takahara, F., Ohira, Y., & Iwasaki, K. 2011, MNRAS, 415, 3434, doi: 10.1111/j.1365-2966.2011.18980.x
  * Gabici et al. (2009) Gabici, S., Aharonian, F. A., & Casanova, S. 2009, MNRAS, 396, 1629, doi: 10.1111/j.1365-2966.2009.14832.x
  * Haslam et al. (1980) Haslam, C. G. T., Pauls, T., & Salter, C. J. 1980, A&A, 92, 57
  * Jeong et al. (2012) Jeong, I.-G., Byun, D.-Y., Koo, B.-C., et al. 2012, Ap&SS, 342, 389, doi: 10.1007/s10509-012-1196-1
  * Jiang et al. (2007) Jiang, B., Chen, Y., & Wang, Q. D. 2007, ApJ, 670, 1142, doi: 10.1086/522299
  * Kaplan et al. (2004) Kaplan, D. L., Frail, D. A., Gaensler, B. M., et al. 2004, ApJS, 153, 269, doi: 10.1086/421065
  * Lalitha et al. (1984) Lalitha, P., Mantovani, F., Salter, C. J., & Tomasi, P. 1984, A&A, 131, 196
  * Lande et al. (2012) Lande, J., Ackermann, M., Allafort, A., et al. 2012, ApJ, 756, 5, doi: 10.1088/0004-637X/756/1/5
  * Landecker et al. (1999) Landecker, T. L., Routledge, D., Reynolds, S. P., et al. 1999, ApJ, 527, 866, doi: 10.1086/308100
  * Li et al. (2023) Li, Y., Xin, Y., Liu, S., & He, Y. 2023, ApJ, 945, 21, doi: 10.3847/1538-4357/acb81d
  * Lorimer et al. (1998) Lorimer, D. R., Lyne, A. G., & Camilo, F. 1998, A&A, 331, 1002, doi: 10.48550/arXiv.astro-ph/9801096
  * Ohira et al. (2011) Ohira, Y., Murase, K., & Yamazaki, R. 2011, MNRAS, 410, 1577, doi: 10.1111/j.1365-2966.2010.17539.x
  * Ptuskin et al. (2006) Ptuskin, V. S., Moskalenko, I. V., Jones, F. C., Strong, A. W., & Zirakashvili, V. N. 2006, ApJ, 642, 902, doi: 10.1086/501117
  * Ptuskin & Zirakashvili (2005) Ptuskin, V. S., & Zirakashvili, V. N. 2005, A&A, 429, 755, doi: 10.1051/0004-6361:20041517
  * Reid et al. (2016) Reid, M. J., Dame, T. M., Menten, K. M., & Brunthaler, A. 2016, ApJ, 823, 77, doi: 10.3847/0004-637X/823/2/77
  * Reid et al. (2019) Reid, M. J., Menten, K. M., Brunthaler, A., et al. 2019, ApJ, 885, 131, doi: 10.3847/1538-4357/ab4a11
  * Roger & Costain (1976) Roger, R. S., & Costain, C. H. 1976, A&A, 51, 151
  * Sedov (1959) Sedov, L. I. 1959, Similarity and Dimensional Methods in Mechanics
  * Straal & van Leeuwen (2019) Straal, S. M., & van Leeuwen, J. 2019, A&A, 623, A90, doi: 10.1051/0004-6361/201833922
  * Thoudam & Hörandel (2012) Thoudam, S., & Hörandel, J. R. 2012, MNRAS, 419, 624, doi: 10.1111/j.1365-2966.2011.19724.x
  * Uchiyama et al. (2012) Uchiyama, Y., Funk, S., Katagiri, H., et al. 2012, ApJ, 749, L35, doi: 10.1088/2041-8205/749/2/L35
  * Wood et al. (2017) Wood, M., Caputo, R., Charles, E., et al. 2017, in International Cosmic Ray Conference, Vol. 301, 35th International Cosmic Ray Conference (ICRC2017), 824.  https://arxiv.org/abs/1707.09551
  * Zabalza (2015) Zabalza, V. 2015, Proc. of International Cosmic Ray Conference 2015, 922
  * Zhang et al. (2021) Zhang, H.-M., Liu, R.-Y., Su, Y., et al. 2021, ApJ, 923, 106, doi: 10.3847/1538-4357/ac36c6
  * Zhang & Xin (2023) Zhang, P., & Xin, Y. 2023, ApJ, 951, 142, doi: 10.3847/1538-4357/acdaca