# Unifying Local and Global Multimodal Features for Place Recognition
in Aliased and Low-Texture Environments

Alberto García-Hernández1,2, Riccardo Giubilato1, Klaus H. Strobl1, Javier
Civera2, and Rudolph Triebel1 This work was supported by the Helmholtz
Association, project ARCHES (contract no. ZT-0033), the Spanish Government
(projects PID2021-127685NB-I00 and TED2021-131150B-I00) and the Aragón
Government (DGA-T45-23R).1Institute of Robotics and Mechatronics, German
Aerospace Center (DLR<EMAIL_ADDRESS>Universidad de
Zaragoza, Spain<EMAIL_ADDRESS>
###### Abstract

Perceptual aliasing and weak textures pose significant challenges to the task
of place recognition, hindering the performance of Simultaneous Localization
and Mapping (SLAM) systems. This paper presents a novel model, called UMF
(standing for Unifying Local and Global Multimodal Features) that 1) leverages
multi-modality by cross-attention blocks between vision and LiDAR features,
and 2) includes a re-ranking stage that re-orders based on local feature
matching the top-k candidates retrieved using a global representation. Our
experiments, particularly on sequences captured on a planetary-analogous
environment, show that UMF outperforms significantly previous baselines in
those challenging aliased environments. Since our work aims to enhance the
reliability of SLAM in _all_ situations, we also explore its performance on
the widely used RobotCar dataset, for broader applicability. Code and models
are available at https://github.com/DLR-RM/UMF.

##  I Introduction

Simultaneous Localization and Mapping (SLAM) has emerged as a central
technology in a multitude of industries including autonomous driving [1, 2],
automated construction [3], and agriculture [4, 5]. Its development and
adoption have been accelerated by advancements in sensor technologies,
including multi-camera setups, RGB-D sensors, and more recently, 3D Light
Detection and Ranging (LiDAR) sensors, facilitating the construction of large-
scale, dense, high-resolution, and consistent 3D maps.

Many SLAM benchmarks, such as KITTI [6], Oxford RobotCar [7], KAIST [8], and
4Season [9], are focused on autonomous driving in urban environments. While
these datasets incorporate substantial challenges, the highly structured
scenarios and vehicle-centric perspectives may oversimplify odometry and place
recognition tasks. Other datasets like TUM RGB-D [10] and TUM VI [11] offer
sequences captured using hand-held stereo or RGB-D cameras primarily indoors.
However, these sequences are typically short and re-visit previous locations
from similar viewpoints and with limited variations in visual appearance. To
address these limitations, synthetic datasets such as ICL-NUIM [12] and
TartanAir [13] simulate more general and challenging motions and environments.
However, in the quest for more diversity, unstructured natural environments
still pose the most challenging conditions for visual or LiDAR-based SLAM.

(a) The LRU rover traversing a planetary-like environment.

(b) Aligned visual and 3D LiDAR data.

Figure 1: (a) The LRU rover traversing the Moon-analogue environment of Mt.
Etna, Sicily, recording the DLR Planetary Stereo Solid-State LiDAR Inertial
(S3LI) dataset [14]. (b) Aligned visual and 3D LiDAR data. Note the
challenging texture and geometry for place recognition.

In this paper we propose a novel multimodal place recognition method that we
denote as UMF, standing for Unifying local and global Multimodal
Features.111Our naming builds on the title of Ref. [15]. Our model leverages
deep local and global features from visual and LiDAR data, fusing both
modalities via cross-attention mechanisms. As our main novelty, we incorporate
local feature-based re-ranking to a multimodal setup, showing in our
experiments that it leads to a substantial and consistent improvement, in
particular in challenging unstructured natural scenes.

##  II Related Work

Place recognition is a well-established field with a wide literature, relevant
applications, and still significant challenges [16, 17, 18, 19, 20, 21, 22].
In relation to our research, a first categorization can be done based on the
input data modality. Visual place recognition has largely evolved from
traditional hand-crafted features [23, 24] to learned descriptors [25, 26],
recently moving towards attention mechanisms [27, 28]. LiDAR place recognition
emphasizes the geometric properties of environments, generally utilizing point
cloud data. PointNet [29, 26] serves as the foundational architecture for
subsequent works such as PointNetVLAD [30] and LPD-Net [11]. Attention
mechanisms have also been used to improve feature specificity [31, 32].
Visual-LiDAR place recognition aims at overcoming the limitations inherent in
using either modality alone, combining their strengths for more reliable
recognition. Techniques like feature concatenation and attention are used in
[33, 34], respectively, which are the main baselines in this area. Attention
mechanisms and transformers have shown great utility in several ways. They can
serve as patch descriptor filters to spotlight important scene elements [35,
27] or as weight maps that modulate the CNN feature maps to generate global
features [36]. This enables models to focus on salient features and ignore
irrelevant information.

Patch-NetVLAD [37] and the work of Zheng et al. [38] showcase the benefits of
unifying local and global features for visual place recognition. The latter
targets the problem of visual aliasing in unstructured scenes [39], combining
SuperPoint and SuperGlue in the context of 6D camera relocalization. In this
work, we propose for the first time unifying local and global features for
multimodal data.

Self-supervised pretraining refers to an initial training stage across
multiple domains that can then be fine-tuned for specific downstream tasks.
Here, contrastive learning methods like SimCLR [40] and BYOL [41] take on an
important role. Likewise, generative self-supervised methods allow to learn
the underlying dependencies in data. In particular, the masked autoencoders do
so by learning to fill artificially corrupted parts of the input data, which
can only be successfully performed by producing rich, useful data
representations. Unlike contrastive methods, generative methods do not require
carefully designed data augmentations or pair constructions, which are
especially challenging on point clouds. For visual data, ConvNeXt [42] and
Spark [43] employ 3D sparse convolutions to apply CNNs to point clouds. For
LiDAR data, however, the application of self-supervised learning presents
unique challenges. The inherent sparsity and lack of order in LiDAR point
clouds, along with the necessity to capture complex geometric and spatial
features, make conventional generative methods less suitable. One of the few
successful works for generative, self-supervised pretraining on large point
clouds is Occupancy-MAE [44]. It employs voxelization and masking during the
training phase, prompting the model to become “voxel-aware,” thereby
efficiently leveraging the geometrical and spatial redundancies within the
point cloud data. Occupancy-MAE has proven effective in downstream tasks,
including 3D object detection and semantic segmentation, even with a high
masking ratio of up to 70%. In this work, we employ generative self-supervised
pretraining methods in the spirit of the masked autoencoder due to their
suitability for multiple data modalities and their lack of dependence on
intricate 3D data augmentations. In particular, we build from the Occupancy-
MAE due to its proven capability to manage large-scale point clouds with a
high masking ratio.

##  III Unifying Local and Global
Multimodal Features (UMF)

Figure 2: UMF overview. Each branch encodes each of the inputs independently.
The encodings of the individual modalities are fused by self- and cross-
attention modules into a single global multimodal representation. For each
individual data modality, separate branches extract also local features.
During inference, we query a database of places with the global multimodal
descriptor using a K-Dimensional Tree, the top-$k$ candidates are retrieved
via NN-search, and finally they are re-ranked using local features from both
modalities. This last stage is the main contribution of our paper.

An overview of our model is showcased in Fig. 2. At its core, the UMF design
consists of 1) two distinct branches, one for visual data and the other for
LiDAR, that encode the data and extract local features, 2) self- and cross-
attention blocks to extract a global representation leveraging both
modalities, and 3) re-ranking methods that take the top-$k$ candidates by
their global representation and re-order them based on local feature matching
on each modality. Next, we provide further details on the main aspects of the
method.

Image and LiDAR Encoding: Firstly, each of these branches transform the
respective input data into low-dimensional representations via a Resnet50
backbone and a LiDAR convolutional encoder [44]. Both encoders follow a
Feature Pyramid Network (FPN) architecture to extract multi-scale features,
thereby capturing both the local fine-grained details and also relevant
patterns with wider extent.

Self- and Cross-Attention: Following [34], the first work that incorporated
attention in multimodal data for place recognition, our UMF model also
incorporates attention to enhance its capability to dynamically focus on
different parts of the input data. In the self-attention layers, the model
assigns different weights of importance to the features within a single
modality (either $F_{\mathrm{Vision}}$ or $F_{\mathrm{LiDAR}}$), thereby
capturing patterns within local and global contexts. This allows the model to
identify distinctive patterns within each modality and enhances its ability to
recognize places based on a single modality. Cross-attention layers, on the
other hand, take features from both modalities ($F_{\mathrm{Vision}}$ and
$F_{\mathrm{LiDAR}}$) as inputs. By interleaving self- and cross-attention
layers within our UMF model, it becomes capable of capturing relevant patterns
between the two modalities, thereby learning richer scene representations.

Local and Global Features: Building upon [15], our UMF model incorporates both
local and global features, integrating fine-grained details along with the
global spatial embedding, and extends the approach to the multi-modal case. We
utilize transformers with positional encoding for coarse-level fusion,
employing both self- and cross-attention. We implement the same approach
introduced in DETR [45], which ensures that each element in the feature maps
$F_{\mathrm{Vision}}$ and $F_{\mathrm{LiDAR}}$ has a unique positional
encoding such that the transformed features exhibit dependencies on relative
positions. The relative position dependency enhances the model spatial
awareness, as well as inter- and intra-modality relationships between features
in the fusion branch. Furthermore, such a design consideration promotes
rotation and viewpoint invariance, enhancing robustness against perspective
changes.

As the main novelty of our work, we incorporate re-ranking strategies to
multimodal place recognition models. Specifically, we evaluate two strategies
for re-ranking based on matching local features, the first one using the so-
called Super-features [28] and the second one implementing RANSAC geometric
verification.

Super-Features: Super-features [28] were proposed as mid-level scene
representations, showing excellent results in place recognition tasks. The
Local Super-features Integration Transformer (LIT) is trained via contrastive
learning, passing local features through a transformer layer, as summarized in
Fig. 3. Specifically, pairs of Super-features are constrained by the
contrastive margin loss $\mathcal{L}_{s}$, which minimizes the pairwise
distance between matching pairs while simultaneously reducing the spatial
redundancy of Super-features $s$ within an image:

$\mathcal{L}_{s}=\sum_{\mathscr{P}}\left[\left\|s-s^{+}\right\|_{2}^{2}+\sum_{n}\left[\mu^{\prime}-\|s-n\|_{2}^{2}\right]^{+}\right],$
(1)

where $\mu$ is a margin hyper-parameter. The index $n$ corresponds to the
Super-features extracted from all negative images in the training tuple, which
are compared against a specific Super-feature $s$. More explicitly, $n$ serves
as a collection of negative samples with the same Super-feature ID as that of
$s$ (denoted as $i(s)$).

This process results in a $N\\!\times\\!F$ ordered set of $N$ Super-features
of $F$ dimensions. The construction of Super-features involves an iterative
attention module, generating a set where each element focuses on a localized
and discriminative image pattern. In this paper, we extend the approach of
[28] to 3D scenarios, extracting Super-features from images as well as from
voxelized point clouds.

Figure 3: Local Super-features extracted with the LIT module for both
modalities. Attention maps show the areas where each Super-feature is focused
on.

To create Super-features that are complementary, they are encouraged to attend
to different ($i\\!\neq\\!j$) local features, i.e., to different image
locations. To do that, the cosine similarity between the attention maps of all
Super-features of every image is minimized. Let matrix
$\boldsymbol{\alpha}=\left[\tilde{\boldsymbol{\alpha}}_{1},\ldots,\tilde{\boldsymbol{\alpha}}_{N}\right]$
denote the $N$ attention maps after the last iteration of LIT. The attention
decorrelation loss is then given by:

$\mathcal{L}_{\mathrm{attn}}(\boldsymbol{x})=\frac{1}{N(N-1)}\sum_{i\neq
j}\frac{\tilde{\boldsymbol{\alpha}}_{i}^{\top}\cdot\tilde{\boldsymbol{\alpha}}_{j}}{\left\|\tilde{\boldsymbol{\alpha}}_{i}\right\|_{2}\left\|\tilde{\boldsymbol{\alpha}}_{j}\right\|_{2}}\quad.$
(2)

In our UMF model, we implement a re-ranking mechanism to leverage
correspondences at Super-feature level. As mentioned above, for any Super-
feature $s\\!\in\\!\mathcal{S}$ we have a function $i(s)$ that returns the
Super-feature ID, i.e., $i\left(s_{i}\right)\\!=\\!i,\forall
s_{i}\\!\in\\!\mathcal{S}$. Also, let
$n(s,\delta)\\!=\\!\arg\min_{s_{i}\in\delta}\left|s-s_{i}\right|_{2}$ be the
nearest neighbor of $s$ from the set $\delta$.

Given a positive pair of images $\boldsymbol{x},\boldsymbol{x}^{+}$, and two
Super-features $\\{s\\!\in\\!\delta,\,s^{\prime}\\!\in\\!\delta^{\prime}\\}$
from their respective Super-feature sets $\\{\delta,\delta^{\prime}\\}$, we
impose the following criteria to consider the Super-feature pair
$\\{s,s^{\prime}\\}$ eligible:

  1. 1.

Reciprocal nearest neighbors: $s\\!=\\!n(s^{\prime},\mathcal{\delta})$ and
$s^{\prime}\\!=\\!n(s,\mathcal{S}^{\prime})$.

  2. 2.

Pass Lowe’s first-to-second neighbor ratio test [46]:
$\left|\boldsymbol{s}-\boldsymbol{s}^{\prime}\right|_{2}/\left|\boldsymbol{s}^{\prime}-n\left(\boldsymbol{s}^{\prime},\delta\backslash{\boldsymbol{s}}\right)\right|_{2}\geqslant\tau$.

  3. 3.

Have the same Super-feature ID: $i(s)=i\left(s^{\prime}\right)$.

We set $\tau=0.8$ after empirical considerations.

RANSAC: The Random Sample Consensus (RANSAC) variant of our UMF model
emphasizes salient feature selection and geometric verification. Initially, a
transformer layer is applied to generate attention maps. These maps then
filter and pinpoint salient local features, retaining those that carry
significant information for place recognition, as shown in Fig. 4. The local
features are further processed using stacked layers of transformer blocks and
we use the attention maps to filter as salient those features for which
$\mathcal{L}_{\mathrm{attn}}\\!>\\!\delta$, $\delta$ being a hyperparameter
that we tuned experimentally for both modalitites. Similarly to the Super-
features, we project the local features to a 1-D embedding and use a
constrastive loss on each local modality during training.

Figure 4: RANSAC local branch visualization. The resulting attention maps are
used to select the salient features.

The RANSAC algorithm is then employed to estimate the geometric transformation
between the current observation and candidate locations. This provides a
robust approach to matching local features, accommodating the presence of
noise and outliers effectively. The spatial consistency score is given by the
number of inliers returned when fitting a homography between the two images or
computing the rigid transformation between voxel grids, using corresponding
keypoints computed by nearest neighbor matching.

###  III-A Training Pipeline

UMF leverages unlabeled data from similar domains, such as the Mars-analogue
in Morocco [39], for pretraining. This self-supervised learning approach makes
the encoder robust to environmental variations, minimizing the dependency on
labeled data and accelerating convergence during fine-tuning for the
downstream task. In detail, it leverages the masked autoencoders for both
visual and LiDAR modalities, inspired by Spark [43] and Occupancy-MAE [44].

Visual self-supervised pretraining follows a contrastive learning approach on
a vast unlabeled dataset. It is designed to discern visually similar yet
distinct locations, thereby improving the model’s ability to tackle visual
aliasing. A patch-wise masking strategy segments images into non-overlapping
square patches, each subject to independent masking according to a
predetermined mask ratio. We follow the approach of [43], where the authors
propose to assemble all unmasked patches into a sparse image, which is then
encoded using sparse convolutions. The encoder $f$ is based on ResNet and the
decoder $g$ on U-Net, including three blocks with upsampling layers (see Fig.
5). Subsequently, a “densification” process, based on mask embeddings,
produces dense feature and projection layers.

Figure 5: Illustration of our pre-training on RobotCar. First, masked inputs
are encoded by $f$, followed by the densification process in the decoder $g$.
After pre-training, only the encoder $f$ is used for downstream tasks.

LiDAR self-supervised pretraining adopts the self-supervised pre-training
scheme based on masked autoencoders on voxelized clouds presented in Ref.
[47]. This encourages the network to reason about high-level semantics to
recover the masked occupancy distribution of the 3D scene from a limited
number of visible voxels. A binary occupancy classification loss is calculated
using cross-entropy:

$\mathcal{L}_{\mathrm{LiDAR}}=-\frac{1}{n_{\mathrm{batch}}}\sum_{i=1}^{n_{\mathrm{batch}}}\sum_{j=1}^{n_{\mathrm{voxel}}}\mathbf{T}_{j}^{i}\log\mathbf{P}_{j}^{i}\quad,$
(3)

where $\mathbf{P}_{j}^{i}$ is the predicted occupancy probability of voxel $j$
for the $i$-th training sample and $\mathbf{T}_{j}^{i}$ the ground truth
probability indicating whether the voxel contains points.

To address the problem of the range-dependent density of LiDAR points, we
employ a range-aware random masking strategy [44]. This separates the occupied
voxels into three groups based on their distance from the sensor. We apply a
distinct random masking strategy to each group, with decreasing ratio with
increasing distance.

Downstream task fine-tuning for place recognition uses a triplet margin loss
and batch-hard negative mining, following the MinkLoc approach [33]. The
triplets are constructed based on spatial proximity, with a 12 m radius for
similarity and over 60 m for dissimilarity. Batch-hard negative mining targets
active triplets for effective model refinement. Properly balancing global and
modality-specific losses is here of major importance.

##  IV Experimental Results

We test the UMF model on two very different datasets: the DLR S3LI Dataset for
planetary rover exploration [14] and the Oxford RobotCar dataset for
autonomous driving [48].

###  IV-A Data Setup

The Oxford RobotCar dataset include diverse driving scenarios across varied
weather and lighting conditions. Nonground point clouds are down-sampled to
4096 points, while corresponding RGB images are down-sampled from
1280$\times$960 to 224$\times$224\. To enhance data diversity and limit
overfitting, we randomly sample from 15 closest RGB images during training,
while only one RGB image with the closest timestamp is used during evaluation.
Similarity is defined based on their spatial proximity: elements within 10 m
are deemed similar, those separated by at least 50 m are dissimilar; and those
falling between 10 and 50 m are neutral. The dataset is split into disjoint
training (21.7k elements) and test (3k elements) areas based on UTM
coordinates, following the evaluation protocol and the train/test split
(baseline) introduced in [33].

The DLR S3LI dataset includes sequences captured with a hand-held sensor setup
comprising a solid-state LiDAR and a stereo camera. The planetary environment
is affected by extreme visual aliasing and lacks of salient visual or
structural features. The dataset was split into training and validation set,
using s3li_loops and s3li_traverse_1 for testing, as they contain overlapping
areas.

Additional datasets for pre-training and fine-tuning  were used to increase
the model’s robustness and performance. For unstructured planetary
environments, we incorporate the MADMAX dataset [39], the Erfoud dataset [49],
and the Long Range Navigation Tests (LRNTs) [50]. To overcome the scarcity of
datasets that contain both visual and LiDAR modalities, we also generated
synthetic sequences with the OAISYS photorealistic simulator [51].

Dataset | # Samples | Modalities used
---|---|---
MADMAX [39] | 11,000 | Stereo
Erfoud [49] | 7,000 | Stereo, LiDAR
LRNTs [50] | 5,600 | Stereo
OAISYS [51] | 10,000 | Stereo, LiDAR, Instance, Semantic
TABLE I: Datasets used for pre-training and fine-tuning.

###  IV-B Implementation Details

All models were implemented in PyTorch [52] and trained on a compute cluster
equipped with 8 NVIDIA RTX 3090 GPUs. For fine-tuning, the learning rate was
set to $1\mathrm{e}\\!-\\!5$ and reduced by a factor of $1\mathrm{e}\\!-\\!1$
upon plateauing. The model was trained for 200 epochs using Adam [53]. The
input image size was set to 224$\times$224\. When re-ranking global feature
retrieval results with local feature-based matching, the top 25 ranked images
from the first stage are considered.

Super-features: The score for Super-features was determined using
$\mathcal{L}_{\mathrm{attn}}$. The attention maps were generated at image size
of 56$\times$56 and voxel size of 50$\times$50$\times$50\. Super-features were
represented as a tensor of size $[N,F]$, and the dimensions for the visual and
point cloud features are set to 128 and 32, respectively. The final ranking is
based on the number of matching features that satisfy the criteria described
in Sec. III.

RANSAC: We use a multi head transformer encoder in our model to process the
fine features. The model returns the average of all attention maps and selects
an optimal threshold to identify keypoints. The output consists of attention
maps for the image and voxel, of sizes $[N,56,56]$ and $[N,50,50,50]$,
respectively, where $N$ is the number of keypoints. The feature maps for the
image and voxel are $[56,56,128]$ and $[50,50,50,32]$, respectively. The re-
ranking of candidates is based on the total number of inliers, either for one
or both modalities:
$score=\\#\mathrm{inliers}_{\mathrm{Vision}}+\\#\mathrm{inliers}_{\mathrm{LiDAR}}$.

Similarity Threshold: The threshold $\theta$ plays a pivotal role in our
model’s ranking process, acting as a cut-off value to distinguish between
similar and dissimilar pairs during re-ranking. Given two samples $S_{1}$ and
$S_{2}$, with a derived similarity score $\mathrm{sim}(S_{1},S_{2})$, they are
classified as recognizing the same place if
$\mathrm{sim}(S_{1},S_{2})\\!>\\!\theta$, where $\theta$ adjusts the model
sensitivity when classifying candidate matches as true or false positives, and
it is determined empirically through cross-validation to optimize precision
and recall rates.

###  IV-C Comparison against Baselines

Table II shows the comparative performance of our UMF models vs. relevant
baselines on the S3LI dataset. Our best model outperforms the best baseline by
more than 2% in the three metrics chosen, which can be attributed mainly to
the use of local features for re-ranking. While the LiDAR data has a limited
field of view, they still provide a valuable input under challenging
conditions by reducing uncertainty. Moreover, the accurate depths play a
crucial role in establishing correct positive pairs. Also, the geometric
verification using RANSAC shows a significant impact in aliased environments,
outperforming other approaches. For an illustrative understanding, Fig. 6
depicts the Super-features attention maps generated by our LIT.

Method | Recall@1 | Recall@5 | Top 1% recall
---|---|---|---
Visual
DBoW2 | 37.44 | 66.1 | 68.12
NetVLAD | 67.2 | 75.5 | 78.3
MinkLoc++ (Vision) | 68.8 | 77.3 | 79.2
LiDAR
PointNet++ | 48.41 | 67.77 | 71.8
MinkLoc++ (LiDAR) | 42.4 | 65.8 | 69.4
Multimodal
MinkLoc++ | 71.4 | 80.1 | 85.2
AdaFusion | 73.1 | 82.3 | 87.2
UMF (only global feat.) | 73.5 | 82.9 | 87.5
UMF (Super-features) | 75 | 85.1 | 89.1
UMF (RANSAC) | 75.3 | 85.3 | 89.5
TABLE II: Comparison against baselines on S3LI (Mt. Etna), for single- and
multi-modality. Figure 6: Illustration of the Super-features attention maps
generated by the Learning Iterative Transformer (LIT), for S3LI (1st rows) and
RobotCar (2nd and 3rd row). The first three Super-features are highlighted,
showing the model’s propensity to focus on semantic patterns such as rock
formations and terrain structures. The last row shows the attention map
extracted from the multi head attention for a pair of RobotCar samples.

The evaluation on the RobotCar dataset, in a very different application domain
and with different sensor specifications, demonstrates the robustness and
adaptability of UMF. Table III shows multimodal results compared against the
baselines MinkLoc++ and AdaFusion, in which our UMF model with re-ranking
outperforms again both baselines.

Method | Recall@1 | Top 1% recall
---|---|---
Multimodal
MinkLoc++ | 96.7 | 99.1
AdaFusion | 98.1 | 99.2
UMF (only global feat.) | 97.9 | 99.1
UMF (Super-features) | 98.1 | 99.1
UMF (RANSAC) | 98.3 | 99.3
TABLE III: Comparison against baselines on RobotCar. | S3LI (Mt. Etna) | RobotCar
---|---|---
Modality | Method | Recall@1 | Recall@5 | Top 1% recall | Recall@1 | Recall@5 | Top 1% recall
Vision | UMF (Super-features) | 74.5 | 84.3 | 89 | 98 | 98.4 | 99.1
UMF (RANSAC) | 75.1 | 84.9 | 89.3 | 98.1 | 98.5 | 99.2
LiDAR | UMF (Super-features) | 73.7 | 83.4 | 87.5 | 97.8. | 98.2 | 99.1
UMF (RANSAC) | 73.9 | 83.8 | 87.8 | 98 | 98.3 | 99.1
Vision+LiDAR | UMF (Super-features) | 75 | 85.1 | 89.1 | 98.1 | 98.5 | 99.1
UMF (RANSAC) | 75.3 | 85.3 | 89.5 | 98.3 | 98.8 | 99.3
TABLE IV: Ablation study on S3LI and RobotCar. Note how the multimodal
(Vision+LiDAR) setup presents the best metrics, effectively leveraging both
modalities. Observe also how our RANSAC variant consistently outperforms
Super-features. Finally, note how the improvement offered by our methods is
bigger in S3LI than in the almost saturated RobotCar.

###  IV-D Ablation Studies

Re-ranking: To analyze the role of the re-ranking module in our approach, we
conduct experiments involving variations of it, such as with and without, for
one or all modalities. The improvements after re-ranking shown in Table IV
substantiate its importance. On the S3LI dataset, we observe that the visual
modality outperforms LiDAR. The re-ranking step on LiDAR data shows marginal
improvements. The fusion of both modalities, however, does aid in overcoming
visual challenges such as poor lighting or aliasing in both datasets. Notably,
RANSAC emerged as a clear winner over Super-features within the re-ranking
methods. This advantage, however, comes at the expense of an increased
computational cost.

The normalized similarity threshold $\theta$ defined in Section IV-B has also
been adjusted to examine the effectiveness of each variant in the precision-
recall curves, see Fig. 7. As expected, RANSAC clearly outperforms the Super-
features curve. Despite both approaches offering competitive performance
compared to the baseline models without re-ranking, precision deteriorates
rapidly, confirming again challenges in S3LI due to the lack of salient
features.

Figure 7: Precision-recall curves on S3LI (Mt. Etna). We compare the base UMF
model with both re-ranking variants.

Pre-training: To assess the impact of pre-training (see Sec. III-A), we
perform an ablation study by initializing the weights with either random
values or a pre-trained set. Table V shows the important role of pre-training.
In our experience, this especially applies for complex models. Also, it
reinforces the model’s robustness and generalization by reducing its
propensity to overfit on the visual modality.

Method | Recall@1 | Recall@5 | Top 1% recall
---|---|---|---
UMF (w/o pre-training) | 70.4 | 81.2 | 85.9
UMF (w/ pre-training) | 73.5 | 82.9 | 87.5
TABLE V: Influence of pre-training on S3LI.

##  V Conclusions

In this paper we have presented the “Unifying local and global Multimodal
Features” (UMF) model, a novel place recognition method that fuses local and
global features of both visual and LiDAR data using transformers and
incorporates re-ranking steps based on single-modality local features. We
evaluate our UMF against state-of-the-art baselines on two different domains:
urban and planetary. UMF shows superior performance in terms of Recall@$\\!N$
in both domains. In particular, we observed that our model outperforms
significantly previous baselines in the extreme conditions of the S3LI dataset
(Mt. Etna, Sicily), demonstrating the potential of a sound multimodal fusion
for place recognition in challenging (unstructured and aliased) scenes.

## References

  * [1] G. Bresson _et al._ , “Simultaneous localization and mapping: A survey of current trends in autonomous driving,” _IEEE Transactions on Intelligent Vehicles_ , vol. 2, no. 3, pp. 194–220, 2017.
  * [2] M. Burki _et al._ , “Vizard: Reliable visual localization for autonomous vehicles in urban outdoor environments,” in _IEEE Intelligent Vehicles Symposium (IV)_ , 2019, pp. 1124–1130.
  * [3] S. Mascaro _et al._ , “Towards automating construction tasks: Large-scale object mapping, segmentation, and manipulation,” _Journal of Field Robotics_ , vol. 38, no. 5, pp. 684–699, 2021.
  * [4] X. Shu _et al._ , “Slam in the field: An evaluation of monocular mapping and localization on challenging dynamic agricultural environment,” in _IEEE/CVF Winter Conference on Applications of Computer Vision_ , 2021, pp. 1761–1771.
  * [5] R. Oliveira _et al._ , “Advances in agriculture robotics: A state-of-the-art review and challenges ahead,” _Robotics_ , vol. 10, no. 2, p. 52, 2021.
  * [6] N. Gelfand _et al._ , “Geometrically stable sampling for the icp algorithm,” in _International Conference on 3-D Digital Imaging and Modeling (3DIM)_ , 2003, pp. 260–267.
  * [7] A. Geiger _et al._ , “Vision meets robotics: The kitti dataset,” _The International Journal of Robotics Research_ , vol. 32, no. 11, pp. 1231–1237, 2013.
  * [8] J. Choi _et al._ , “Kaist multi-spectral day/night data set for autonomous and assisted driving,” _IEEE Transactions on Intelligent Transportation Systems_ , vol. 19, no. 3, pp. 934–948, 2018.
  * [9] S. Wenzel _et al._ , “4seasons: A cross-season dataset for multi-weather slam in autonomous driving,” in _DAGM German Conference on Pattern Recognition_. Springer, 2020, pp. 404–417.
  * [10] J. Sturm _et al._ , “A benchmark for the evaluation of rgb-d slam systems,” in _IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_ , 2012, pp. 573–580.
  * [11] R. Schubert _et al._ , “The tum vi benchmark for evaluating visual-inertial odometry,” in _IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_ , 2018, pp. 1680–1687.
  * [12] T. Schops _et al._ , “Bad slam: Bundle adjusted direct rgb-d slam,” in _IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR)_ , 2019\.
  * [13] A. Handa _et al._ , “A benchmark for rgb-d visual odometry, 3d reconstruction and slam,” in _IEEE International Conference on Robotics and Automation (ICRA)_ , 2014, pp. 1524–1531.
  * [14] R. Giubilato, W. Stürzl, A. Wedler, and R. Triebel, “Challenges of slam in extremely unstructured environments: The dlr planetary stereo, solid-state lidar, inertial dataset,” _IEEE Robotics and Automation Letters_ , vol. 7, no. 4, pp. 8721–8728, 2022.
  * [15] B. Cao, A. Araujo, and J. Sim, “Unifying deep local and global features for image search,” in _Computer Vision–ECCV 2020: 16th European Conference, Glasgow, UK, August 23–28, 2020, Proceedings, Part XX 16_. Springer, 2020, pp. 726–743.
  * [16] B. Williams, M. Cummins, J. Neira, P. Newman, I. Reid, and J. Tardós, “A comparison of loop closing techniques in monocular slam,” _Robotics and Autonomous Systems_ , vol. 57, no. 12, pp. 1188–1197, 2009.
  * [17] E. Garcia-Fidalgo and A. Ortiz, “Vision-based topological mapping and localization methods: A survey,” _Robotics and Autonomous Systems_ , vol. 64, pp. 1–20, 2015.
  * [18] S. Lowry, N. Sünderhauf, P. Newman, J. J. Leonard, D. Cox, P. Corke, and M. J. Milford, “Visual place recognition: A survey,” _ieee transactions on robotics_ , vol. 32, no. 1, pp. 1–19, 2015.
  * [19] X. Zhang, L. Wang, and Y. Su, “Visual place recognition: A survey from deep learning perspective,” _Pattern Recognition_ , vol. 113, p. 107760, 2021\.
  * [20] C. Masone and B. Caputo, “A survey on deep visual place recognition,” _IEEE Access_ , vol. 9, pp. 19 516–19 547, 2021.
  * [21] S. Garg, T. Fischer, and M. Milford, “Where is your place, visual place recognition?” _arXiv preprint arXiv:2103.06443_ , 2021.
  * [22] S. Schubert, P. Neubert, S. Garg, M. Milford, and T. Fischer, “Visual place recognition: A tutorial,” _arXiv preprint arXiv:2303.03281_ , 2023.
  * [23] D. Gálvez-López and J. D. Tardós, “Bags of binary words for fast place recognition in image sequences,” _IEEE Transactions on Robotics_ , vol. 28, no. 5, pp. 1188–1197, October 2012.
  * [24] C. Campos, R. Elvira, J. J. G. Rodríguez, J. M. M. Montiel, and J. D. Tardós, “ORB-SLAM3: An accurate Open-Source library for visual, Visual–Inertial, and multimap SLAM,” _IEEE Trans. Rob._ , pp. 1–17, 2021.
  * [25] R. Arandjelović, P. Gronat, A. Torii, T. Pajdla, and J. Sivic, “NetVLAD: CNN architecture for weakly supervised place recognition,” in _IEEE Conference on Computer Vision and Pattern Recognition_ , 2016.
  * [26] C. R. Qi, L. Yi, H. Su, and L. J. Guibas, “Pointnet++: Deep hierarchical feature learning on point sets in a metric space,” _arXiv preprint arXiv:1706.02413_ , 2017.
  * [27] Y. Zhu, J. Wang, L. Xie, and L. Zheng, “Attention-based pyramid aggregation network for visual place recognition,” 2018.
  * [28] P. Weinzaepfel, T. Lucas, D. Larlus, and Y. Kalantidis, “Learning Super-Features for image retrieval,” in _International Conference on Learning Representations_ , 2022.
  * [29] C. R. Qi, H. Su, K. Mo, and L. J. Guibas, “Pointnet: Deep learning on point sets for 3d classification and segmentation,” _arXiv preprint arXiv:1612.00593_ , 2016.
  * [30] M. A. Uy and G. H. Lee, “Pointnetvlad: Deep point cloud based retrieval for large-scale place recognition,” in _The IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ , 2018.
  * [31] T.-X. Xu, Y.-C. Guo, Y.-K. Lai, and S.-H. Zhang, “TransLoc3D: Point cloud based large-scale place recognition using adaptive receptive fields,” May 2021\.
  * [32] J. Ma, J. Zhang, J. Xu, R. Ai, W. Gu, and X. Chen, “OverlapTransformer: An efficient and Yaw-Angle-Invariant transformer network for LiDAR-Based place recognition,” _IEEE Robotics and Automation Letters_ , vol. 7, no. 3, pp. 6958–6965, July 2022.
  * [33] J. Komorowski, M. Wysoczańska, and T. Trzcinski, “MinkLoc++: Lidar and monocular image fusion for place recognition,” in _2021 International Joint Conference on Neural Networks (IJCNN)_ , July 2021, pp. 1–8.
  * [34] H. Lai, P. Yin, and S. Scherer, “AdaFusion: Visual-LiDAR fusion with adaptive weights for place recognition,” _IEEE Robotics and Automation Letters_ , vol. 7, no. 4, pp. 12 038–12 045, 2022.
  * [35] M.-H. Guo, C.-Z. Lu, Z.-N. Liu, M.-M. Cheng, and S.-M. Hu, “Visual attention network,” 2022.
  * [36] R. Wang, Y. Shen, W. Zuo, S. Zhou, and N. Zheng, “TransVPR: Transformer-based place recognition with multi-level attention aggregation,” Jan. 2022.
  * [37] S. Hausler, S. Garg, M. Xu, M. Milford, and T. Fischer, “Patch-NetVLAD: Multi-Scale fusion of Locally-Global descriptors for place recognition,” 2021\.
  * [38] Y. Zheng, T. Birdal, F. Xia, Y. Yang, Y. Duan, and L. J. Guibas, “6d camera relocalization in visually ambiguous extreme environments,” _arXiv preprint arXiv:2207.06333_ , 2022.
  * [39] L. Meyer, M. Smíšek, A. Fontan Villacampa, L. Oliva Maza, D. Medina, M. J. Schuster, F. Steidle, M. Vayugundla, M. G. Müller, B. Rebele, _et al._ , “The madmax data set for visual-inertial rover navigation on mars,” _Journal of Field Robotics_ , vol. 38, no. 6, pp. 833–853, 2021.
  * [40] T. Chen, S. Kornblith, K. Swersky, M. Norouzi, and G. Hinton, “Big self-supervised models are strong semi-supervised learners,” _arXiv preprint arXiv:2006.10029_ , 2020.
  * [41] J.-B. Grill, F. Strub, F. Altché, C. Tallec, P. H. Richemond, E. Buchatskaya, C. Doersch, B. A. Pires, Z. D. Guo, M. G. Azar, B. Piot, K. Kavukcuoglu, R. Munos, and M. Valko, “Bootstrap your own latent: A new approach to self-supervised learning,” 2020.
  * [42] S. Woo, S. Debnath, R. Hu, X. Chen, Z. Liu, I. S. Kweon, and S. Xie, “ConvNeXt v2: Co-designing and scaling ConvNets with masked autoencoders,” Jan. 2023.
  * [43] K. Tian, Y. Jiang, Q. Diao, C. Lin, L. Wang, and Z. Yuan, “Designing BERT for convolutional networks: Sparse and hierarchical masked modeling,” Jan. 2023\.
  * [44] C. Min, X. Xu, D. Zhao, L. Xiao, Y. Nie, and B. Dai, “Occupancy-mae: Self-supervised pre-training large-scale lidar point clouds with masked occupancy autoencoders,” _arXiv e-prints_ , 2022.
  * [45] N. Carion, F. Massa, G. Synnaeve, N. Usunier, A. Kirillov, and S. Zagoruyko, “End-to-end object detection with transformers,” in _Computer Vision–ECCV 2020: 16th European Conference, Glasgow, UK, August 23–28, 2020, Proceedings, Part I 16_. Springer, 2020, pp. 213–229.
  * [46] D. G. Lowe, “Distinctive image features from scale-invariant keypoints,” _International journal of computer vision_ , vol. 60, pp. 91–110, 2004.
  * [47] G. Hess, J. Jaxing, E. Svensson, D. Hagerman, C. Petersson, and L. Svensson, “Masked autoencoder for self-supervised pre-training on lidar point clouds,” in _2023 IEEE/CVF Winter Conference on Applications of Computer Vision Workshops (WACVW)_. IEEE, jan 2023. [Online]. Available: https://doi.org/10.1109
  * [48] D. Barnes, M. Gadd, P. Murcutt, P. Newman, and I. Posner, “The oxford radar robotcar dataset: A radar extension to the oxford robotcar dataset,” in _2020 IEEE International Conference on Robotics and Automation (ICRA)_. IEEE, 2020, pp. 6433–6438.
  * [49] S. Lacroix, A. De Maio, Q. Labourey, E. P. Mendes, P. Narvor, V. Bissonette, C. Bazerque, F. Souvannavong, R. Viards, and M. Azkarate, “The erfoud dataset: a comprehensive multi-camera and lidar data collection for planetary exploration,” in _15th Symposium on Advanced Space Technologies in Robotics and Automation_ , 2020.
  * [50] M. Vayugundla, F. Steidle, M. Smisek, M. J. Schuster, K. Bussmann, and A. Wedler, “Datasets of long range navigation experiments in a moon analogue environment on mount etna,” in _ISR 2018; 50th International Symposium on Robotics_ , 2018, pp. 1–7.
  * [51] M. G. Müller, M. Durner, A. Gawel, W. Stürzl, R. Triebel, and R. Siegwart, “A Photorealistic Terrain Simulation Pipeline for Unstructured Outdoor Environments,” in _IEEE/RSJ International Conference on Intelligent Robots and Systems_ , 2021.
  * [52] F. A. research lab, “Pytorch,” _https://pytorch.org/_ , (date accessed 11-3-2023).
  * [53] D. P. Kingma and J. Ba, “Adam: A method for stochastic optimization,” _arXiv preprint arXiv:1412.6980_ , 2014.


# Broadband X-ray timing and spectral characteristics of the accretion-powered
millisecond X-ray pulsar MAXI J1816$-$195

Zhaosheng Li Key Laboratory of Stars and Interstellar Medium, Xiangtan
University, Xiangtan 411105, Hunan, China  Lucien Kuiper SRON-Netherlands
Institute for Space Research, Niels Bohrweg 4, 2333 CA, Leiden, The
Netherlands  Mingyu Ge Key Laboratory of Particle Astrophysics, Institute of
High Energy Physics, Chinese Academy of Sciences, 19B Yuquan Road, Beijing
100049, China  Maurizio Falanga International Space Science Institute (ISSI),
Hallerstrasse 6, 3012 Bern, Switzerland  Physikalisches Institut, University
of Bern, Sidlerstrasse 5, 3012 Bern, Switzerland  Juri Poutanen Department of
Physics and Astronomy, University of Turku, FI-20014, Finland  Long Ji School
of Physics and Astronomy, Sun Yat-sen University, Zhuhai, 519082, People’s
Republic of China  Yuanyue Pan  Key Laboratory of Stars and Interstellar
Medium, Xiangtan University, Xiangtan 411105, Hunan, P.R. China  Yue Huang
Key Laboratory of Particle Astrophysics, Institute of High Energy Physics,
Chinese Academy of Sciences, 19B Yuquan Road, Beijing 100049, China  Renxin Xu
Department of Astronomy, School of Physics, Peking University, Beijing 100871,
China  Kavli Institute for Astronomy and Astrophysics, Peking University,
Beijing 100871, China  Liming Song Key Laboratory of Particle Astrophysics,
Institute of High Energy Physics, Chinese Academy of Sciences, 19B Yuquan
Road, Beijing 100049, China  Jinlu Qu  Key Laboratory of Particle
Astrophysics, Institute of High Energy Physics, Chinese Academy of Sciences,
19B Yuquan Road, Beijing 100049, China  Shu Zhang  Key Laboratory of Particle
Astrophysics, Institute of High Energy Physics, Chinese Academy of Sciences,
19B Yuquan Road, Beijing 100049, China  Fangjun Lu  Key Laboratory of Particle
Astrophysics, Institute of High Energy Physics, Chinese Academy of Sciences,
19B Yuquan Road, Beijing 100049, China  Shuang-Nan Zhang Key Laboratory of
Particle Astrophysics, Institute of High Energy Physics, Chinese Academy of
Sciences, 19B Yuquan Road, Beijing 100049, China

###### Abstract

We studied the broadband X-ray timing and spectral behaviors of the newly
confirmed accreting millisecond X-ray pulsar MAXI J1816$-$195 during its 2022
outburst. We used the data from Insight-HXMT ME/HE, NICER and NuSTAR which
cover the energy range between 0.8–210 keV. A coherent timing analysis of
solely Insight-HXMT HE data across the full outburst revealed a complex
behavior of the timing residuals, also prominently visible in independent
Insight-HXMT ME and NICER data, particularly at rising part of the outburst
and at the very end in NICER data. Therefore, we broke down the full outburst
into a (noisy) rising part, covering only about five days from MJD 59737.0 to
59741.9, and a decaying part lasting for 19 days across MJD 59741.9–59760.6.
Fitting for the decaying part a timing model including a frequency $\nu$ and
frequency time derivative $\dot{\nu}$ component yielded a value of $(+9.0\pm
2.1)\times 10^{-14}~{}{\rm Hz~{}s^{-1}}$ for $\dot{\nu}$, which could be
interpreted as a spin-up under our model assumptions. We detected the X-ray
pulsations up to $\sim$95 keV in a combination of Insight-HXMT HE
observations. The pulse profiles were quite stable over the whole outburst and
could be well described by a truncated Fourier series using two harmonics, the
fundamental and the first overtone. Both components kept alignment in the
range 0.8–64 keV. The joint and time-averaged NICER and Insight-HXMT spectra
in the energy range 1–150 keV were well fitted by the absorbed Comptonization
model compps plus disk blackbody with two additional Gaussian components.
Using the bolometric flux and spin-up values both evaluated during the decay
phase, we determined a magnetic field strength of $(0.2-2)\times 10^{8}$ G for
MAXI J1816$-$195.

pulsars: individual: MAXI J1816$-$195 – stars: neutron – X-rays: general –
X-rays: binaries

††facilities: Insight-HXMT, NICER, NuSTAR

##  1 Introduction

As a subclass of neutron star low-mass X-ray binary (NS LMXB), accreting
millisecond X-ray pulsars (AMXPs) are usually confirmed by its coherent
pulsation during outburst, with spin periods of a few milliseconds and orbital
periods between 40 min and 5 h based on the current sample (see Campana & Di
Salvo, 2018; Di Salvo & Sanna, 2022; Papitto et al., 2020; Patruno & Watts,
2021, for reviews). The recycling scenario suggested that AMXPs are
progenitors of (radio) millisecond pulsars (Alpar et al., 1982; Radhakrishnan
& Srinivasan, 1982). This picture has been verified by the discoveries of
AMXPs starting with SAX J1808.4–3658 by RXTE (Wijnands & van der Klis, 1998),
discovery of the pulsar spin-up induced by accretion (Falanga et al., 2005a),
and by swinging between accretion- and rotation-powered states in IGR
J18245–2452 (Papitto et al., 2013).

The new X-ray transient MAXI J1816$-$195 was discovered by MAXI/Gas Slit
Camera on June 7 during its 2022 outburst (Negoro et al., 2022). Using the
Swift localization, NICER detected pulsations and thermonuclear type I X-ray
bursts from MAXI J1816$-$195, confirming that the source is an AMXP with a
spin frequency of 528 Hz, an orbital period of 4.83 hr and a projected semi-
major axis of 0.26 lt-s (Kennea et al., 2022; Bult et al., 2022b, c, a). The
mass of the companion star was determined to be in the range 0.10–0.55
$M_{\odot}$ (Bult et al., 2022a). Meanwhile, Insight-HXMT found X-ray
pulsations from MAXI J1816$-$195 in the hard X-ray/soft $\gamma$-ray band (Li
et al., 2022). Bult et al. (2022a) proposed the flux-bias model, which
accounts for the effects of the accretion torque and/or the wandering of the
hot spot on the NS surface, to explain the structures left in the timing
residuals across the outburst. Chen et al. (2022) reported the detection of 73
type I X-ray bursts by Insight-HXMT from MAXI J1816$-$195, and obtained an
upper limit of the distance as 6.3 kpc. Wang et al. (2023, in preparation)
found the type I X-ray bursts rate decreasing for the bolometric flux higher
than $\sim 1.04\times 10^{-8}\rm{erg~{}cm^{-2}~{}s^{-1}}$ and estimated the
distance to the source at $\sim 3.4$ kpc (Cavecchi et al., 2020).

Since the onset of the 2022 outburst a search for radio, optical and near-
infrared emission from MAXI J1816$-$195 had been carried out. The radio
emission of MAXI J1816$-$195 showed significant evolution, at 1.28 GHz with a
flux density of $4.2\pm 0.4~{}{\rm mJy}$ on 2022 June 9 and a $3\sigma$ upper
limit of $70{~{}\rm\mu Jy}$ on 2022 June 11 by MeerKAT, at 5.5 GHz with a flux
density of $135\pm 35~{}{\rm\mu Jy}$ on 2022 June 12 by the Australia
Telescope Compact Array (ATCA), and at 6 GHz with a flux density of $457\pm
20~{}{\rm\mu Jy}$ on 2022 June 18 by VLA (Beauchamp et al., 2022; Bright et
al., 2022). However, radio pulsations were not detected. The likely optical
and near-infrared counterpart of MAXI J1816$-$195 had been identified (de
Martino et al., 2022; Kong, 2022). When MAXI J1816$-$195 evolved towards the
quiescent state about one month after the start of the outburst, X-ray
emission from a dust scattered halo around the source position was detected by
Swift (Beardmore et al., 2022).

In this work, we present the broad band timing and spectral analysis of MAXI
J1816$-$195\. We introduce the instruments, Insight-HXMT, NICER and NuSTAR,
and the performed observations in Sect. 2. The timing and the broad-band
spectral results are reported in Sect. 3 and Sect. 4, respectively, and our
findings are discussed in Sect. 5.

##  2 Observations and data reduction

We analyzed the data from MAXI J1816$-$195 collected by Insight-HXMT, NICER,
and NuSTAR.

###  2.1 Insight-HXMT

Insight-HXMT (Insight Hard X-ray Modulation Telescope, Zhang et al., 2020) is
the first Chinese X-ray telescope, and is equipped with three slat-collimated
instruments: the Low Energy X-ray telescope (LE, 1–12 keV; Chen et al., 2020),
the Medium Energy X-ray telescope (ME, 5–35 keV; Cao et al., 2020) and the
High Energy X-ray telescope (HE, 20–350 keV; Liu et al., 2020), providing
capabilities for broadband X-ray timing and spectroscopy. Insight-HXMT carried
out high-cadence observations of MAXI J1816$-$195 starting on MJD 59738.17,
$\sim$0.5 days after the detection by MAXI. The set of 79 ME and 77 HE
observations includes runs P0404275001 – P0404275024. The HE and ME data were
studied to investigate the timing and spectral properties of the source
because of their broadband X-ray coverage and good time resolution ( $\sim$2
$\mu$s for HE, $\sim$20 $\mu$s for ME). We analyzed the data using the
Insight-HXMT Data Analysis Software (HXMTDAS) version 2.05. The ME and HE data
were calibrated by using the scripts mepical and hepical, respectively. The
good time intervals were individually selected from the scripts megtigen and
hegtigen with the standard criteria, that is, ELV $>$ 0, the satellite located
outside the South Atlantic Anomaly region and the offset angle from the
pointing direction is smaller than $0\farcs 04$. Totally, 73 type I X-ray
bursts have been observed by ME (Chen et al., 2022), which were removed in the
timing and spectral analysis. No bursts were found in HE because the burst
spectra are usually soft. Background subtracted light curves for the ME and HE
were generated (see Figure 2.3). The outburst light curve of Insight-HXMT ME
showed a fast rise during the early days of the outburst and a slow decay in
the next 25 d. However, the Insight-HXMT HE light curve for the 20–50 keV band
dropped continuously from $\sim 100$ cts s-1 since the start of the
observations to the quiescent level during the outburst. The spectra and their
response matrix files are produced by the tools hespecgen and herspgen for HE,
and mespecgen and merspgen for ME, respectively. Finally, we obtained the
cleaned events using mescreen and hescreen and barycentered with the tool
hxbary.

###  2.2 NICER

NICER started regular observations of MAXI J1816$-$195 on 2022 June 7 (MJD
59737.6; obs. ID 5202820101), several hours after the detection by MAXI, and
ended these on 2022 August 3 (MJD 59794.01; obs. ID 5533013504) when the
source was already several weeks in its quiescent state. We performed the
standard data processing using the NICER Data Analysis Software (NICERDAS).
The default filtering criteria were applied to extract the cleaned event data.
Next, we extracted 1-s light curves by xselect for the 0.5–10 and 12–15 keV
bands to search for type I X-ray bursts and flaring background, respectively.
Time intervals containing X-ray bursts or background flares were removed in
further analysis. The 0.5–10 keV light curve showed a profile similar to
Insight-HXMT ME, i.e., a fast rise in first few days, and a decay to quiescent
level during the next $\sim 25$ days (see the top panel of Figure 2.3).

The background spectra were produced from the nibackgen3C50 tool (Remillard et
al., 2022). The redistribution matrix file and ancillary response file were
generated from the tools nicerarf and nicerrmf, respectively.

In the pulse-profile analysis of the NICER data (see Sect. 3.2) we used obs.
IDs 5533013001–5533013502 (MJD 59769.566–59786.993; cleaned exposure time
14.388 ks) performed when the source had already reached its quiescent state
with undetectable (pulsed) flux levels, for the estimation of the local
background to obtain background corrected fractional amplitudes.

###  2.3 NuSTAR

NuSTAR observed MAXI J1816$-$195 on 2022 June 23 (obs. ID 90801315001; MJD
$\sim 59753.456-59754.413$) for about 35.7 ks. We cleaned the event file using
the NuSTAR pipeline tool nupipeline for both FPMA and FPMB. The light curves
and spectra were extracted from a circular region with a radius of 100″
centered on the source location by using nuproducts, and the response and
ancillary response files were produced simultaneously. To extract the
background spectra, we chose a source-free circular background region centered
at
$(\alpha_{2000},\delta_{2000})=(18^{\rm{\scriptsize{h}}}16^{\rm{\scriptsize{m}}}28\fs
4048,-19\arcdeg 39\arcmin 56\farcs 308)$ with a radial aperture of
$100\arcsec$. From the light curves, four type I X-ray bursts were identified
(see also Mandal et al., 2023). For the timing and spectral analysis the time
intervals during these bursts are discarded resulting in a total exposure time
of 29.6 ks.

In the timing analysis, we barycentered event data from the source region,
adopting a $90\arcsec$ extraction radius, using HEASOFT multi-mission tool
barycorr v2.16, with NuSTAR fine clock-correction file #142, yielding time
tags accurate at 60–100 $\mu$s level (Bachetti et al., 2021) in absolute time.
To obtain background corrected timing characteristics, such as fractional
amplitudes, we used the above mentioned background region, however, now with a
90″ extraction radius.

The light curves of MAXI J1816$-$195 during its 2022 outburst. From top to
bottom: the 100 s binned light curves from NICER, Insight-HXMT ME, and
Insight-HXMT HE are displayed, respectively. The energy range of each light
curve is indicated in each panel. The grey area marks the time of NuSTAR
observation. Time intervals covering X-ray bursts are removed. The blue dashed
lines represent the time interval that we used for the derivation of our
timing models, see Table 1 and Sect. 3.

##  3 Timing analysis

Irrespective the instrument in timing analyses we converted the Terrestial
Time (TT) arrival times of the selected events to arrival times at the solar
system barycenter (in a Barycentric Dynamical Time (TDB) time scale). For this
process, throughout in this work we used: (1) the JPL DE405 solar system
ephemeris; (2) the instantaneous spacecraft position with respect to the
Earth’s center; (3) the (most accurate) VLA position of MAXI J1816$-$195,
$\alpha_{\rm 2000}=18^{\rm h}16^{\rm m}52\fs 41168$ and $\delta_{\rm
2000}=-19\arcdeg 37\arcmin 57\farcs 40138$ (Beauchamp et al., 2022), which
differ by $0\farcs 4$ in declination from the Swift location (Kennea et al.,
2022). Moreover, in the selection procedures we ignored time intervals during
which type I X-ray bursts or high-background flaring occurred.

###  3.1 Insight-HXMT and NICER Timing analysis

Initially, we constrained the timing analysis to Insight-HXMT HE data using
only events with measured energies between 20 and 35 keV. Before generating
phase-coherent timing models (ephemerides) describing accurately the spin
frequency evolution across the outburst we first need to correct the
barycentered event arrival times for the periodic orbital motion effects.
Early mission data covering only a part (MJD 59738–59752) of the full dataset
available for Insight-HXMT HE analysis had been used to improve the orbital
parameters of MAXI J1816$-$195 applying an optimisation scheme based on a
SIMPLEX algorithm (see De Falco et al., 2017a; Li et al., 2021, for more
details) by iteratively improving the $Z_{2}^{2}(\phi)$-test statistics
(Buccheri et al., 1983) with respect to the pulse frequency $\nu$, the orbital
period $P_{\rm orb}$, and time of the ascending node $T_{\rm asc}$. The
optimised values for the (most critical) orbital parameters $P_{\rm orb}$ and
$T_{\rm asc}$, while fixing the eccentricity to zero and the projected semi-
major axis $a_{\rm x}\sin i$ to the value quoted in Bult et al. (2022b) are
listed in Table 1. The obtained values are fully consistent within $1\sigma$
with the best orbital parameters obtained later by Bult et al. (2022a) using
the full NICER database. Therefore, we kept the values for the orbital
parameters as derived from our three parameter optimisation scheme using 20–35
keV Insight-HXMT HE data.

Table 1: Positional, orbital and spin parameters derived in this work and used as fixed values from literature for MAXI J1816$-$195. Parameter | Values | Units
---|---|---
$\alpha_{2000}$ | $18^{\hbox{\scriptsize h}}16^{\hbox{\scriptsize m}}52\fs 41168(12)$ |
$\delta_{2000}$ | $-19\arcdeg 37\arcmin 57\farcs 40138(473)$ |
JPL Ephemeris | DE405 |
$P_{\rm orb}$ | $17402.5786(60)$ | s
$a_{\rm x}\,\sin i$ | $0.262\,948(18)$ | lt-s
$e$ |  $0$ (fixed) |
$T_{\rm asc}$ | $59738.875\,632(4)$ | MJD (TDB)
Constant Frequency model
Validity range | $59741.9-59760.6$ | MJD (TDB)
$t_{0}$ (Epoch) | $59741.0$ | MJD (TDB)
$\nu$ | $528.611\,105\,832(4)$ | Hz
$\chi^{2}$/d.o.f | $86.78/(47-1)=1.886$ |
| Quadratic / Spin-up model |
Validity range | $59741.9-59760.6$ | MJD (TDB)
$t_{0}$ (Epoch) | $59741.0$ | MJD (TDB)
$\nu$ | $528.611\,105\,774(12)$ | Hz
$\dot{\nu}$ | $(9.0\pm 2.1)\times 10^{-14}$ | $\rm{Hz~{}s^{-1}}$
$\chi^{2}$/d.o.f | $68.18/(47-2)=1.515$ |

After correcting for the orbital motion effects we performed a time-of-arrival
(ToA; see Kuiper & Hermsen, 2009, for more details) analysis to obtain
accurate estimates for the pulse frequency $\nu$, and (if required) its first
time derivative $\dot{\nu}$ over a certain time interval. We obtained 69 ToA’s
for Insight-HXMT HE (20–35 keV) covering MJD 59738–59761. Similar analyses
yielded 73 ToA’s for Insight-HXMT ME restricted to the 10–20 keV band (MJD
59738–59761), and 26 ToA’s for NICER using 1.8–10 keV events across MJD
59737–59761.

Figure 1: Pulse phase residuals of Insight-HXMT HE (open black triangles;
20–35 keV), Insight-HXMT ME (open orange squares; 10–20 keV) and NICER (open
aqua diamonds; 1.8–10 keV) ToAs with respect to a two-parameter timing model
(including a quadratic / spin-up term) valid for the range MJD
59741.9–59760.6. These models are based on solely Insight-HXMT HE ToAs (black
data points). The interval bounded by the two blue dashed lines represents
this time interval that has been used to derive the timing models given in
Table 1. Note, that there are no ToA measurements beyond MJD 59760.6 for
Insight-HXMT HE (and ME) contrary to NICER.

Fitting the complete Insight-HXMT HE ToA dataset with a constant frequency
model resulted in a very poor fit with an unacceptable $\chi^{2}$. Ignoring 22
ToA’s from the first $\sim 5$ days of the outburst yielded reasonable fits
assuming a constant frequency model and a model including a frequency
derivative term. Both are listed in Table 1. The constant frequency model is
fully consistent at $1\sigma$ confidence with the corresponding model, based
on solely NICER data, shown in Table 1 of Bult et al. (2022a) for 46 degrees
of freedom.

The model with both the $\nu$ and $\dot{\nu}$ terms, however, improves upon
the constant frequency model by reducing the $\chi^{2}$ from 86.78 to 68.18
adding one fit parameter. Applying a maximum likelihood ratio test to assess
the model improvement yields an improvement of
$(86.78-68.18)^{0.5}=4.3\sigma$, adopting one degree of freedom (=the
additional fit parameter).111The square-root relation holds only for the case
of one additional fit parameter (= one degree of freedom) The measured
$\dot{\nu}$ of $(9.0\pm 2.1)\times 10^{-14}~{}\rm{Hz~{}s^{-1}}$ can be
interpreted as a spin-up across the time interval MJD 59741.9–59760.6. We have
also investigated the effects of employing a flux-bias model in the ToA fit
procedure by adding a (flux) bias term in the form,
$\Phi_{\rm{\scriptsize{bias}}}=\alpha(C/C_{\max})^{\Gamma}$ with
$C_{\max}=895.37\ {\rm c\,s^{-1}}$ (the maximum observed count rate in the
0.5–10 keV NICER bandpass), to the pulse-phase residual relation in analogy
with Bult et al. (2022a) (see their Eq. 1). Note, that in Table 1 of Bult et
al. (2022a) both optimized flux-bias model parameters, their $b$ and $\Gamma$,
are listed abusively with a plus sign. Taken this sign inconsistency into
account we could reproduce their findings and obtained $\Gamma=-1.05\pm 0.15$,
consistent within one sigma confidence with their value of $-1.20\pm 0.20$
employing the full NICER ToA dataset running from MJD 59737.748–59762.610.
Constraining the NICER ToA dataset to MJD 59742.370–59762.610, and so ignoring
the rising and maximum part of the light curve (see Figure 2.3), we obtained a
flux-bias index of $\Gamma=-0.77\pm 0.09$, however, as with the full range
NICER model shown before, with a statistically unacceptable fit quality of
$\chi^{2}_{r}\sim 3.54$ compatible with the results shown in Bult et al.
(2022a). Next, we fitted a constant frequency- plus flux-bias model with a
fixed index $\Gamma=-0.77$ to the Insight-HXMT-HE ToA dataset covering MJD
59741.9–59760.6, the decaying part of the outburst. This model improves upon
the constant frequency model by reducing the $\chi^{2}$ from 86.78 to 83.10
adding one free parameter (the flux-bias model scale $\alpha$), which
represents a $(86.78-83.10)^{0.5}=1.92\sigma$ improvement with respect to the
constant frequency model, applying a maximum likelihood ratio-test. The best
fit flux-bias scale parameter is $\alpha=(-7.1\pm 3.7)\times 10^{-3}$, while
$\nu=528.611\,105\,851(11)$ Hz. Therefore, the spin-up model for the decaying
part of the outburst (covering MJD 59741.9–59760.6) provides a better data
description than the flux-bias model for the Insight-HXMT-HE ToA dataset.

The pulse-phase residuals of the complete sets of ToA’s from all three
instruments combined (each determined for a different energy band) after
folding upon the (Insight-HXMT HE 20–35 keV based) spin-up model are shown in
Figure 1. Irrespective the energy band each instrument showed a quite similar
behaviour: during the first $\sim 5$ days of the outburst - the rising flux
part -there is a structured residual feature which can be described by a
(short-duration) spin-down model. This feature is the cause of the very poor
$\chi^{2}$ when the full ToA datasets are considered (see also Table 1 of Bult
et al., 2022a, indicating unacceptable $\chi^{2}$ fit values). Also, the steep
positive excursions for the NICER ToAs at the early beginning and very end
(not covered by Insight-HXMT observations) of the outburst are striking in
Figure 1. It is clear that even a (bolometric) flux bias model as well as the
constant frequency and quadratic models considered in this work cannot
describe the full set of residuals, consistent with Bult et al. (2022a).

Figure 2: The 0.8–210 keV broadband pulse-phase distributions of MAXI
J1816$-$195 observed by NICER (panels a–c, 0.8–10 keV), NuSTAR (panels g–k,
3–79 keV), and Insight-HXMT (panels d–f, 5–35 keV for ME; panels l–p, 20–210
keV for HE). In order to improve clarity, two cycles are shown. The error bars
represent $1\sigma$ errors. The morphology is almost unchanged with energy.
All profiles reach their maximum near phase $\sim$0.5.

Figure 3: Left: The fractional amplitude of various harmonics: 1 – fundamental
(black crosses) and 2 – first overtone (red crosses) as a function of energy
using background-corrected data from NICER (filled squares), NuSTAR
(triangles), Insight-HXMT ME (filled circles) and Insight-HXMT HE (open
circles). Right: The phase angle (in units radians/$2\pi$) as a function of
energy for the two Fourier components, the fundamental (in black) and the
first overtone (in red), which are independent of the background subtraction
method.

###  3.2 Pulse profile analysis

With the accurate orbital and (spin-up) timing models (see Table 1) we
converted the event arrival times from each high-energy instrument to pulse
phases. To obtain high-statistics pulse-phase distributions (pulse profiles
including background) we combined the phase information for each instrument
across the validity interval (MJD 59741.9–59760.9) and sorted these on
(measured) energy. We verified for each instrument that the pulse shape for
each individual observation constituting the combination remained the same
within statistical uncertainties. The measured pulse-phase distributions for
the four different high-energy instruments are shown in Figure 2 for various
different energy bands covering the 0.8–210 keV range.

Panels (a)–(c) of Figure 2 show the NICER pulse profiles for the energy
intervals 0.8–1.8, 1.8–3, and 3–10 keV, respectively, while panels (d)–(f)
show the Insight-HXMT ME profiles for the 5–10, 10–20 and 20–35 keV bands. The
NuSTAR profiles are displayed in panels (g)–(k) covering the 3–79 keV energy
range. Noteworthy is that comparing the NICER and NuSTAR 3–10 keV
distributions, shown in panels (c) and (g), respectively, a small phase shift
of $0.065\pm 0.020$ is detectable which is just compatible with the absolute
timing accuracy of NuSTAR (Bachetti et al., 2021). Pulsed emission is detected
significantly up to $\sim$60 keV, with a non-uniformity significance of
$5.4\sigma$ for the 35–60 keV band (panel j) adopting a $Z_{2}^{2}$-test
(Buccheri et al., 1983). This test applied to the 60–79 keV band (panel k)
yielded an insignificant result of $1.3\sigma$. Finally, panels (l)–(p) show
the Insight-HXMT HE profiles for energies of 20–30, 30–50, 50–65, 65–95 and
95–210 keV, respectively. Interestingly, pulsed emission is detected
significantly up to $\sim$95 keV. The 65–95 keV band (panel o) shows a non-
uniform significance of $4.4\sigma$ applying a $Z_{2}^{2}$-test with
indications for a second weaker pulse near phase 0.15 that seems to pop-up at
energies above $\sim$50 keV (panels n–o). The non-uniformity significance for
the distribution shown in panel (p) for the 95–210 keV band is about
$1.1\sigma$, and so insignificant.

To obtain quantitative information about morphology changes of the pulse
profiles as a function of energy, we fit the pulse profiles with a truncated
Fourier series given by a formula

$F(\phi)=A_{0}+\sum_{k=1}^{2}A_{k}\ \cos[2\pi\ k(\phi-\phi_{k})],$ (1)

where $A_{1}$ and $A_{2}$ are the amplitudes, $\phi_{1}$ and $\phi_{2}$ are
the phase angles (in radians/$2\pi$), of the fundamental and the first
overtone, respectively, and $A_{0}$ is the constant level of the profile. We
verified that higher harmonics are not required to accurately describe the
data. The left panel of Figure 3 shows the background corrected fractional
amplitudes $A_{k}/A_{0}$ of the fundamental (black) and first overtone (red)
for all four instruments covering the 0.8–95 keV band, across which we
detected significant pulsed emission.

The background correction determination was straightforward for NICER and
NuSTAR. We used for the former non-imaging instrument the data from
observations when the source was completely Off (see Sect. 2.2) with a
combined Good Time Interval (GTI) exposure of 14.388 ks. Along with the On-
time GTI exposure of 51.735 ks this yielded a scale factor
$s_{\scriptsize{\rm{bg}}}$ of $3.5957$ to be applied in the background
subtraction. For the latter (imaging) instrument the data from a different
source-free location in the field-of-view with equal area as source region
($s_{\scriptsize{\rm{bg}}}\equiv 1.0$) was used.

For the non-imaging Insight-HXMT ME/HE (with GTI exposures of 186.171 ks and
132.093 ks, respectively, for the validity interval during the On source
period) the last observations of MAXI J1816$-$195 performed on MJD 59765 (2022
July 5) totaling exposure times of 13.474 ks and 9.022 ks for ME and HE,
respectively, were used to estimate the underlying background, when the source
did not show pulsed emission anymore, but was not completely Off (see, e.g.,
Beardmore et al., 2022, who reported the detection of the source at a
$4.1\sigma$ level in 5.2 ks Swift XRT data taken a day later on MJD 59766).

Thus, the ME/HE event data of MJD 59765 still contain a tiny contribution from
MAXI J1816$-$195, that would yield an over-subtraction when applied
uncorrected in the background subtraction procedure. Demanding that the
background subtraction did not yield a negative value for energy bands at the
high end of the sensitivity window of both instruments (ME $\gtrsim 30$ keV;
HE $\gtrsim 95$ keV) we could estimate the genuine background contribution in
both samples. For the ME/HE, we found that $\sim$9% and $\sim$5% of the events
from the background sample was originating from the ‘almost’ Off MAXI
J1816$-$195, respectively. Taking these fractions into account in the
background subtraction process showed that the $A_{1}/A_{0}$ (and also
$A_{2}/A_{0}$) values obtained for ME/HE matched with NuSTAR values obtained
in the overlapping energy range.

It is clear from the left panel of Figure 3 that the pulsed fraction of the
fundamental (black) gradually increases with energy from $\sim$0.5% to 10%–12%
going from $\sim$0.8 to $\sim$10 keV, where a plateau is reached. Above
$\sim$35 keV a decrease sets in to values of 3%–5% above 50 keV. The first
overtone (red) saturates at $\gtrsim$10 keV to a $\sim$2% value, without a
clear decrease beyond. This means that the $A_{2}/A_{1}$ ratio increases the
higher the energy, implying that the first overtone contributes more
significantly at higher energies.

The right panel of Figure 3 shows the phases of the maxima of the fundamental
(black) and first overtone (red) for the four instruments, specifying the
alignment of the components constituting the pulse profile. It is noteworthy
that these quantities, $\phi_{1}$ and $\phi_{2}$, are independent on the used
background method.

For NuSTAR, which provided a snapshot during the full outburst, we have
subtracted 0.065 (in phase units) -still within its absolute timing accuracy -
from the measured $\phi_{1}$ and $\phi_{2}$ values for alignment purposes with
NICER. It is clear that irrespective the instrument both quantities are stable
within statistical uncertainties till $\sim$50 keV, beyond which the first
overtone moves closer to the fundamental. Combined with the increase of the
$A_{2}/A_{1}$ ratio at those energies this means that a second bump/maximum in
the pulse profile in front of the main maximum becomes more prominent at
$\gtrsim$50 keV energies.

##  4 Spectral analysis

The spectral analysis was carried out using xspec version 12.12.1 (Arnaud,
1996). We first performed the joint spectral evolution analysis for NICER and
Insight-HXMT observations (Sect. 4.1). And then, we studied the quasi-
simultaneous broadband spectra of NICER, NuSTAR and Insight-HXMT (Sect. 4.2).
All spectra are grouped by the tool ftgrouppha. For NICER and NuSTAR spectra,
we applied the optimal binning method (Kaastra & Bleeker, 2016). The Insight-
HXMT spectra were grouped with a minimal signal-to-noise ratio of 3 related to
their backgrounds. All uncertainties of the spectral parameters are provided
at a $1\sigma$ confidence level for a single parameter.

###  4.1 The spectral evolution of joint NICER and Insight-HXMT observations

There were many quasi-simultaneous NICER and Insight-HXMT observations of MAXI
J1816$-$195, which allows us to perform a detailed broadband spectral analysis
for the whole outburst. We found that each NICER observation in the MJD
59738–59761 range can be (partially) covered by 1-6 Insight-HXMT observations.
We combined the Insight-HXMT spectra, and then carried out joint NICER and
Insight-HXMT spectral fitting. We adopted the energy ranges of 1–10 keV for
NICER spectra, 10–35 keV for Insight-HXMT ME and 25–150 keV for Insight-HXMT
HE. We fit all the spectra by using the combination of the thermal
Comptonization model, compps, in the slab geometry (Poutanen & Svensson,
1996), the disk blackbody, diskbb, modified by the interstellar absorption
described by the model, tbabs. This model has been used to fit the broadband
spectra of many AXMPs (see, e.g., Gierliński & Poutanen, 2005; Falanga et al.,
2005a, b, 2011, 2012; De Falco et al., 2017b, a; Kuiper et al., 2020; Li et
al., 2018, 2021). To account for cross-calibration uncertainties between
different instruments a multiplication factor is included in the fits. The
factor is fixed at unity for NICER and set free for Insight-HXMT ME/HE. The
main parameters are the Thomson optical depth across the slab, $\tau_{\rm T}$,
the electron temperature, $kT_{\rm e}$, the temperature of the soft seed
photons, $kT_{\rm bb}$, the normalization factor for the seed black body
photons, $K_{\rm bb}=(R_{\rm km}/d_{10})^{2}$ (with $d_{10}$ being distance in
units of 10 kpc), and the inclination angle $\theta$ (fixed at 60°) between
the slab normal and the line of sight to the observer for the compps model;
the disk temperature, $kT_{\rm diskbb}$ and normalization for the multi-black
body model, and the hydrogen column density, $N_{\rm H}$. Most of NICER
spectra showed emission lines around $\sim$1.6 and $\sim$6.5 keV, therefore,
we added two Gaussian components to account for these. The full model is tbabs
(gaussian+gaussian+diskbb+compps) in xspec. The $\sim$1.6 keV component could
have an instrumental origin, while the $\sim$6.5 keV Gaussian component may
come from the iron K$\alpha$ line emission.

The fitted parameters are shown in Figure 4. Most of the spectra can be well
fitted with reduced $\chi^{2}<1.3$. For two observations, the reduced
$\chi^{2}$ are larger than 1.5, which show no features in the NICER residuals
and possibly caused by the spectral evolution during the long exposures. We
calculated the unabsorbed bolometric flux in the 1–250 keV range by using the
tool cflux. During the outburst, the disk blackbody temperature increased from
0.4 to 0.8 keV in the first few days and then dropped back to $\sim$0.4 keV.
The hydrogen column density did not change much, and the mean value is
$(2.31\pm 0.07)\times 10^{22}~{}{\rm cm^{-2}}$. The optical depth was in the
range 1.7–3.0, which explains the visibility of hard X-ray emissions during
the whole outburst. The electron temperature was around 10–25 keV. The
bolometric flux reached a peak value of $1.02\times 10^{-8}~{}{\rm
erg~{}cm^{-2}~{}s^{-1}}$ during the first few days, and in the next 20 days,
it decayed to $8.5\times 10^{-10}~{}{\rm erg~{}cm^{-2}~{}s^{-1}}$ at MJD
59760.9.

We also replaced the compps component by a cutoff power-law model. This model
provided compatible results compared with the compps model, i.e.,
$N_{\rm{H}}$, the temperature and normalization of disk black body component,
the $\chi^{2}$, and the unabsorbed flux.

Figure 4: The fitted results of joint NICER and Insight-HXMT spectra by the
model constant$\times$tbabs(gaussian + gaussian + diskbb + compps). From top
to bottom, the bolometric flux (in units of $10^{-8}$ erg cm-2 s-1), the
hydrogen column density (in $10^{22}$ cm-2), the electron temperature (in
keV), the optical depth, the temperature of disk blackbody (in keV), and the
reduced $\chi^{2}$ are shown.

###  4.2 Broadband spectra including NuSTAR observation

We noted that the only NuSTAR observation was carried out quasi-simultaneously
with some of NICER and Insight-HXMT observations. We fitted the spectra,
including the Insight-HXMT ME and HE data between MJD 59753.46–59754.40 (obs.
ID P0404275015) with a total exposure of 20.99 ks (ME) and 17.51 ks (HE),
NuSTAR FPMA (29.55 ks), FPMB (29.70 ks) and NICER data starting on MJD
59753.30 (5.21 ks, obs. ID 5533011501). We kept the default energy range of
3–79 keV for NuSTAR FPMA/FPMB, and for the other instruments the same ranges
as in Sect. 4.1. The multiplication factors were also included to account for
the cross-calibration uncertainties between the instruments. The factor was
fixed at unity for the NuSTAR FPMA, and set free for the other instruments. We
first fitted the spectra by using an absorbed cutoff power-law model in
combination with a disk black body model and two additional Gaussians. The
best fitted model yielded $N_{\rm H}\approx 2.30\times 10^{22}$ cm-2, a photon
index $\Gamma\approx 1.54$, a cutoff energy $E_{\rm cut}\approx 24.1$ keV, the
Gaussian at 5.44 keV with a width of 1.32 keV, a disk temperature $kT_{\rm
diskbb}\approx 0.61$ keV, and a reduced $\chi^{2}\sim 1.15$ for a total d.o.f.
of 2300. We then replaced the phenomenological cutoff power-law model by the
thermal Comptonization model compps, which is tbabs
(gaussian+gaussian+diskbb+compps). We obtained a best-fit reduced $\chi^{2}$
of 1.17. We noticed that a reflection feature has been reported by Chauhan et
al. (2022), therefore the reflection parameter, $R$, of compps was set free.
The best-fit reduced $\chi^{2}$ decreased to 1.11 for a total d.o.f. of 2298,
better than the cut-off power-law model. The best-fit parameters are listed in
Table 2 and the broad-band spectrum is shown in Figure 5. The NICER spectrum
contributed the most significant residuals. Considering such a large data set,
this model provides an acceptable fit to the spectra. We obtained a
normalization factor for the seed photons, $K_{\rm bb}=65.9\pm 6.8$, which
implies a size of blackbody emitting region $R_{\rm bb}=0.34\sqrt{K_{\rm
bb}}\approx 3~{}{\rm km}$ for a distance of 3.4 kpc (Wang et al. submitted,
see also Chen et al., 2022). This region is smaller than the whole NS surface
with a typical radius of 12 km. However, the true emitting size should be
larger since part of photons were reflected by the accretion disk. It is
worthy to note that it is difficult to constrain the refection contribution by
only using the joint NICER and Insight-HXMT/ME/HE spectra, i.e., without
NuSTAR data, the $\chi^{2}_{r}$ is 1.0, see Sect. 4.1.

We also changed the compps model to a relativistic reflection model, relxillCp
(Dauser et al., 2016), in which the reflection spectrum is calculated using
the Comptonization continuum nthcomp. After many attempts, during the fit, we
set the iron abundance, $A_{\rm Fe}$, the density of the accretion disk, $\log
N$, the ionization parameter, $\log\xi$, the power-law index of the incident
spectrum, $\Gamma$, the reflection fraction, $R_{f}$, the outer emissivity
index, $q_{2}$, free to vary. The spin parameter was fixed at 0.248 for the NS
spin frequency of 528 Hz. The outer disk radius was fixed at 1000 $R_{\rm
ISCO}$. The inclination angle, $i$, was also fixed at 60 degree. Other
parameters were fixed at their default values. This model gave a best-fit
reduced $\chi^{2}$ of 1.11 for a total d.o.f. of 2295, slightly larger than
the compps model. The results are listed in Table 2.

Table 2:  The spectral fitting to the broadband spectra of MAXI J1816$-$195 with the NICER/NuSTAR/Insight-HXMT observations. The best-fitted results in the middle and right columns are from the model constant$\times$tbabs (gaussian+gaussian+diskbb+compps) and constant$\times$ tbabs (gaussian+gaussian+diskbb+relxillCp), respectively. Parameter (units) | Best-fit values
---|---
$N_{\rm H}~{}(10^{22}{\rm cm}^{-2})$ | $2.25\pm 0.01$ | $2.41\pm 0.01$
$E_{\rm line,~{}1}$ (keV) | $1.62\pm 0.02$ | $1.66\pm 0.02$
$\sigma_{1}$ (keV) | $0.13\pm 0.01$ | $0.09\pm 0.01$
Norm1($\times 10^{-3}$) | $6.5\pm 0.6$ | $4.2\pm 0.7$
$E_{\rm line,~{}2}$ (keV) | $5.53\pm 0.06$ | $5.78\pm 0.05$
$\sigma_{2}$ (keV) | $0.48\pm 0.07$ | $0.76\pm 0.06$
Norm2$(\times 10^{-4})$ | $8\pm 2$ | $26\pm 4$
Inclination (deg) | 60 (fixed) | 60 (fixed)
$kT_{\rm diskbb}$ (keV) | $0.63\pm 0.01$ | $0.55\pm 0.01$
$K_{\rm diskbb}$ (km2) | $770\pm 38$ | $826\pm 58$
$kT_{\rm e}$ (keV) | $35.6\pm 1.1$ | $16.8\pm 0.6$
$q_{2}$ | - | $1.87\pm 0.49$
$\Gamma$ | - | $1.88\pm 0.01$
$\log\xi$ | - | $4.31\pm 0.1$
$\log N$ | - | $17.98\pm 0.38$
$kT_{\rm bb,~{}seed}$ (keV) | $1.31\pm 0.03$ | -
$\tau_{\rm T}$ | $0.89\pm 0.04$ | -
$R$ | $1.34\pm 0.09$ | $6.51\pm 4.67$
Norm${}_{\rm refl.}~{}(\times 10^{-4})$ | - | $8.3\pm 1.2$
$K_{\rm compps}$ (km2) | $65.0\pm 6.8$ | -
$C_{\rm{NuSTAR}/FPMA}$ | 1 (fixed) | 1 (fixed)
$C_{\rm{NuSTAR}/FPMB}$ | $1.01\pm 0.01$ | $1.01\pm 0.01$
$C_{\rm{NICER}}$ | $1.05\pm 0.01$ | $1.05\pm 0.01$
$C_{\rm{Insight-HXMT}/ME}$ | $1.11\pm 0.01$ | $1.11\pm 0.01$
$C_{\rm{Insight-HXMT}/HE}$ | $1.20\pm 0.03$ | $1.21\pm 0.05$
$\chi^{2}/{\rm d.o.f.}$ | 2543.6/2298 | 2550.9/2295
$F_{\rm bol}$ ($10^{-9}$ erg s-1 cm-2)aafootnotemark:  |  $3.97\pm$0.01 |  $3.96\pm$0.01
Figure 5: The unabsorbed spectra of MAXI J1816$-$195 from NICER, NuSTAR FPMA,
NuSTAR FPMB, Insight-HXMT ME and Insight-HXMT HE in the 1–150 keV energy
range. The fitting model is constant$\times$tbabs
(gaussian+gaussian+diskbb+compps). The black dashed lines represent the two
gaussian lines. The best model is plotted as black solid line. The disk
blackbody and compps components are shown in dotted and dash-dot-dot-dot
lines, respectively.

##  5 Summary and discussion

In this paper we reported on detailed timing and spectral analyses of Insight-
HXMT, NICER and NuSTAR data collected from the newly confirmed AMXP MAXI
J1816$-$195 during its 2022 outburst. The outburst of MAXI J1816$-$195 reached
its peak flux within the first few days at $\sim 1\times 10^{-8}~{}{\rm
erg~{}cm^{-2}~{}s^{-1}}$ and decayed into quiescent state in the next $\sim
25$ days. The outburst profile observed from MAXI J1816$-$195 is quite similar
to that of other AMXPs, showing a fast rise ($\sim$1–2 days), an exponential
followed by a linear decay, which can be explained by the disk instability
model (King & Ritter, 1998; Powell et al., 2007).

During the full outburst hard X-ray/soft $\gamma$-ray emission was detected up
to $\sim 200$ keV. The joint NICER and Insight-HXMT spectra in the energy of
1–150 keV were well fitted by a combination of a multi-color disk black body
and Comptonization model compps, with two additional Gaussian lines, exposed
to interstellar extinction. A joint spectral fit including quasi-simultaneous
NuSTAR, NICER and Insight-HXMT ME/HE data showed emission reflected from the
accretion disk, which can be modelled by the reflected compps model or the
relxillCP model.

Usually, the hard X-ray emission from AMXPs is relatively weak and has
typically pulsation amplitudes of less than 10% (see e.g., Papitto et al.
2020, and references therein; see also Falanga et al. 2005a for IGR
J00291+5934 with the pulsed fraction of $\sim 12-20\%$ at $\sim 100$ keV).
Therefore, the pulsations of most AMXPs are mostly detected at soft X-rays
with instruments aboard e.g RXTE, XMM-Newton or NICER. At hard X-rays,
Insight-HXMT observed Swift J1756.9–2508 during its 2018 outburst, and
detected the X-ray pulsations up to $\sim 100$ keV using an ephemeris
determined by NICER at soft X-rays (Li et al., 2021). In the case of MAXI
J1816$-$195, at hard X-rays both the ME and HE aboard of Insight-HXMT were
capable of deriving the pulsar ephemeris directly from their event data,
obtaining parameters well consistent with those derived by NICER at soft
X-rays. Particularly, our timing analysis using the DE405 solar system
ephemeris and the most accurate VLA (Beauchamp et al., 2022) source location
for MAXI J1816$-$195 in the barycentering process, contrary to Bult et al.
(2022a) who used DE430 and the less accurate ($\sim 0\farcs 4$ difference)
Swift location (Kennea et al., 2022), yielded consistent orbital and spin
parameters with Bult et al. (2022a) at $1\sigma$-confidence.

From Insight-HXMT observations hard X-ray pulsations have been detected
significantly up to $\sim 95$ keV. The pulse profiles in 0.8–95 keV band are
well described by a truncated Fourier series. The phase angles of the
fundamental and the first overtone are nearly constant for energies below
$\sim 60$ keV, implying that most photons are likely emitted from the same
region at the NS surface with a very similar emission pattern. The observed
pulse profiles could be used to measure the NS mass and radius (Poutanen &
Gierliński, 2003) of MAXI J1816$-$195 in a future study.

The stability of the pulse profiles across a broad energy range allows us to
determine the absolute timing accuracy of Insight-HXMT in detail. Cross-
correlating Insight-HXMT ME and NICER pulse profiles for the overlapping 5–10
keV energy band yielded a phase shift $\delta\Phi$ of only $-0.008\pm 0.005$,
consistent with zero, suggesting that both instruments are well aligned within
14 $\mu$s, with the Insight-HXMT ME events arriving a little bit later than
NICER ones. A similar correlation analysis now between Insight-HXMT ME and HE
pulse profiles for the 20–35 keV band resulted in a difference of only 0.9
$\mu$s with an uncertainty of $13.9~{}\mu$s.

Bult et al. (2022a) performed a coherent timing analysis adopting a constant
frequency- and a flux bias model using NICER ToA’s covering the full outburst
period. They found unacceptably high $\chi^{2}$ values for both models
assumptions resulting in complex timing residuals for the first 4–5 days
during the rising part of the outburst, a flat shape during the next $\sim 20$
days in the decay phase, and finally, (linearly) increasing phase deviations
at the end of outburst, also visible at the very beginning of the outburst.

In this work we analysed the timing data from Insight-HXMT ME and HE
observations as well as from NICER, and found that the ToAs from the different
(independent) instruments were well consistent with each other in spite their
different bandpass. The timing residuals showed erratic behavior (or V-shaped
wedge) during the rising part of the outburst between MJD 59737.0 and 59741.9,
similar to the results shown in Bult et al. (2022a).

For the decaying part of the outburst between MJD 59741.9–59760.6 we detected
a significant positive quadratic term $\dot{\nu}$ of $(9.0\pm 2.1)\times
10^{-14}~{}\rm{Hz~{}s^{-1}}$ fitting a two-parameter timing model, including
$\nu$ and $\dot{\nu}$, yielding an acceptable $\chi^{2}$, which could
interpreted as a spin-up.

We realize that such spin-up or spin-down terms can arise when an inaccurate
source location is used in the barycentering process. However, in our work we
have used the sub-arcsecond accurate VLA-location for MAXI J1816$-$195, with
an accuracy better than $0\farcs 005$. This uncertainty on the location can
account for an uncertainty in the spin-up rate (see e.g., Sanna et al., 2020,
for the formulae) of only $\sim 2.0\times 10^{-16}~{}{\rm Hz~{}s^{-1}}$, 2.7
order of magnitude smaller than the measured $\dot{\nu}$ value, and much
smaller than our error estimate of $2.1\times 10^{-14}~{}\rm{Hz~{}s^{-1}}$. On
the contrary, if we should have used the Swift-XRT location in the
barycentering process with a centroid $0\farcs 4$ from the VLA-location and an
estimated error radius of $2\farcs 2$ then the uncertainty in the quadratic
term would be $\sim$80 times larger, i.e. of the order of $1.6\times
10^{-14}~{}\rm{Hz~{}s^{-1}}$, comparable to our measurement error. Therefore,
we assume that the measured value for $\dot{\nu}$ during the decay phase is
accurate and genuine, and could be interpreted as a gain of angular momentum
by the NS from the accreted matter.

During the outburst X-ray pulsations have been detected at bolometric flux
states with a flux between $\sim(0.1-1)\times 10^{-8}~{}{\rm
erg~{}cm^{-2}~{}s^{-1}}$. Adopting a source distance of 3.4 kpc, as
constrained by Wang et al. (2023, submitted, see also Chen et al. 2022), these
bolometric fluxes correspond to accretion rates in the range
$(0.012-0.12)\dot{M}_{\rm Edd}$, where $\dot{M}_{\rm Edd}=2\times
10^{-8}~{}M_{\odot}~{}{\rm yr^{-1}}$ is the Eddington rate, assuming that
$L=\eta\dot{M}c^{2}$ and a typical value of $\eta=0.1$ for NS (Frank et al.,
2002). By adopting Equations 11 and 12 and the same assumptions in Psaltis &
Chakrabarty (1999), these limits allow us to constrain the surface dipole
magnetic field strength of $(0.25-74)\times 10^{8}$ G (see also Hartman et
al., 2008). If 20% uncertainty of distance is considered, the magnetic field
is in the range of $(0.2-85)\times 10^{8}$ G. The possible bolometric
correction factor of the flux would not change the magnetic field
significantly because the measured fluxes are covered a broadband energy
range.

During the spin-up period the averaged bolometric flux is $0.56\times
10^{-8}~{}{\rm erg~{}cm^{-2}~{}s^{-1}}$, corresponding to $0.065\dot{M}_{\rm
Edd}$. If the spin-up of the pulsar is solely induced by the angular momentum
transferred from the accreted material to the NS, the magnetic field can be
determined from the relation (Shapiro & Teukolsky, 1983; Tong, 2015; Pan et
al., 2022),

$B=0.1I^{7/2}_{45}R_{6}^{-3}M_{1.4}^{-3/2}\dot{\nu}_{-13}^{7/2}\left(\frac{\dot{M}}{0.065\dot{M}_{\rm
Edd}}\right)^{-3}\times 10^{8}~{}{\rm G},$ (2)

where $I_{45}$, $R_{6}$, $M_{1.4}$, $\dot{\nu}_{-13}$ are the moment of
inertia, the NS radius and mass, the spin up in units of $10^{45}$ g cm2,
$10^{6}$ cm, $1.4M_{\odot}$, and $10^{-13}$ Hz s-1, respectively. By adopting
$I_{45}=1.5$ (see e.g., Worley et al., 2008) and the measured spin-up rate, we
obtain a magnetic field strength of $(0.04-2)\times 10^{8}$ G. Combined with
the aforementioned range, the magnetic field strength of MAXI J1816$-$195 is
$(0.2-2)\times 10^{8}$ G, which is comparable with magnetic field strengths
estimated for several other AMXPs (see e.g., Hartman et al., 2008; Patruno et
al., 2009; Patruno, 2010; Hartman et al., 2009, 2011; Papitto et al., 2011).

##

We appreciate the referee for all valuable comments and suggestions that led
to the improvement of the manuscript. This work was supported the Major
Science and Technology Program of Xinjiang Uygur Autonomous Region (No.
2022A03013-3). Z.S.L., Y.Y.P. and L.J. were supported by National Natural
Science Foundation of China (12103042, U1938107, 12273030, 12173103). S.Z. is
supported by the National Key R&D Program of China (2021YFA0718500). J.P.
acknowledges support from the Academy of Finland grant 333112. This work is
also supported by International Partnership Program of Chinese Academy of
Sciences (Grant No.113111KYSB20190020). This work made use of data from the
Insight-HXMT mission, a project funded by China National Space Administration
(CNSA) and the Chinese Academy of Sciences (CAS), and also from the High
Energy Astrophysics Science Archive Research Center (HEASARC), provided by
NASA’s Goddard Space Flight Center.

## References

  * Alpar et al. (1982) Alpar, M. A., Cheng, A. F., Ruderman, M. A., & Shaham, J. 1982, Nature, 300, 728, doi: 10.1038/300728a0
  * Arnaud (1996) Arnaud, K. A. 1996, in ASP Conf. Ser., Vol. 101, Astronomical Data Analysis Software and Systems V, ed. G. H. Jacoby & J. Barnes, 17
  * Bachetti et al. (2021) Bachetti, M., Markwardt, C. B., Grefenstette, B. W., et al. 2021, ApJ, 908, 184, doi: 10.3847/1538-4357/abd1d6
  * Beardmore et al. (2022) Beardmore, A. P., Kennea, J. A., & Swift/XRT Team. 2022, The Astronomer’s Telegram, 15506, 1
  * Beauchamp et al. (2022) Beauchamp, I., Belvedere, S., Hernandez, M., et al. 2022, The Astronomer’s Telegram, 15481, 1
  * Bright et al. (2022) Bright, J., Russell, T., Tremou, E., et al. 2022, The Astronomer’s Telegram, 15484, 1
  * Buccheri et al. (1983) Buccheri, R., Bennett, K., Bignami, G. F., et al. 1983, A&A, 128, 245
  * Bult et al. (2022a) Bult, P., Altamirano, D., Arzoumanian, Z., et al. 2022a, ApJ, 935, L32, doi: 10.3847/2041-8213/ac87f9
  * Bult et al. (2022b) Bult, P. M., Ng, M., Altamirano, W. I. D., et al. 2022b, The Astronomer’s Telegram, 15425, 1
  * Bult et al. (2022c) Bult, P. M., Sanna, A., Ng, M., et al. 2022c, The Astronomer’s Telegram, 15431, 1
  * Campana & Di Salvo (2018) Campana, S., & Di Salvo, T. 2018, in ASSL, Vol. 457, The Physics and Astrophysics of Neutron Stars, ed. L. Rezzolla, P. Pizzochero, D. I. Jones, N. Rea, & I. Vidaña, 149, doi: 10.1007/978-3-319-97616-7_4
  * Cao et al. (2020) Cao, X., Jiang, W., Meng, B., et al. 2020, Science China Physics, Mechanics, and Astronomy, 63, 249504, doi: 10.1007/s11433-019-1506-1
  * Cavecchi et al. (2020) Cavecchi, Y., Galloway, D. K., Goodwin, A. J., Johnston, Z., & Heger, A. 2020, MNRAS, 499, 2148, doi: 10.1093/mnras/staa2858
  * Chauhan et al. (2022) Chauhan, J., Lohfink, A., Bharali, P., et al. 2022, The Astronomer’s Telegram, 15470, 1
  * Chen et al. (2020) Chen, Y., Cui, W., Li, W., et al. 2020, Science China Physics, Mechanics, and Astronomy, 63, 249505, doi: 10.1007/s11433-019-1469-5
  * Chen et al. (2022) Chen, Y.-P., Zhang, S., Ji, L., et al. 2022, ApJ, 936, L21, doi: https://doi.org/10.3847/2041-8213/ac8c2c
  * Dauser et al. (2016) Dauser, T., García, J., Walton, D. J., et al. 2016, A&A, 590, A76, doi: 10.1051/0004-6361/201628135
  * De Falco et al. (2017a) De Falco, V., Kuiper, L., Bozzo, E., et al. 2017a, A&A, 603, A16, doi: 10.1051/0004-6361/201730600
  * De Falco et al. (2017b) —. 2017b, A&A, 599, A88, doi: 10.1051/0004-6361/201629575
  * de Martino et al. (2022) de Martino, D., D’Avanzo, P., Ambrosino, F., et al. 2022, The Astronomer’s Telegram, 15479, 1
  * Di Salvo & Sanna (2022) Di Salvo, T., & Sanna, A. 2022, in ASSL, Vol. 465, Millisecond Pulsars, ed. S. Bhattacharyya, A. Papitto, & D. Bhattacharya, 87–124, doi: 10.1007/978-3-030-85198-9_4
  * Falanga et al. (2012) Falanga, M., Kuiper, L., Poutanen, J., et al. 2012, A&A, 545, A26, doi: 10.1051/0004-6361/201219582
  * Falanga et al. (2005a) —. 2005a, A&A, 444, 15, doi: 10.1051/0004-6361:20053472
  * Falanga et al. (2005b) Falanga, M., Bonnet-Bidaud, J. M., Poutanen, J., et al. 2005b, A&A, 436, 647, doi: 10.1051/0004-6361:20042575
  * Falanga et al. (2011) Falanga, M., Kuiper, L., Poutanen, J., et al. 2011, A&A, 529, A68, doi: 10.1051/0004-6361/201016240
  * Frank et al. (2002) Frank, J., King, A., & Raine, D. J. 2002, Accretion Power in Astrophysics: Third Edition, 398
  * Gierliński & Poutanen (2005) Gierliński, M., & Poutanen, J. 2005, MNRAS, 359, 1261, doi: 10.1111/j.1365-2966.2005.09004.x
  * Hartman et al. (2011) Hartman, J. M., Galloway, D. K., & Chakrabarty, D. 2011, ApJ, 726, 26, doi: 10.1088/0004-637X/726/1/26
  * Hartman et al. (2009) Hartman, J. M., Patruno, A., Chakrabarty, D., et al. 2009, ApJ, 702, 1673, doi: 10.1088/0004-637X/702/2/1673
  * Hartman et al. (2008) —. 2008, ApJ, 675, 1468, doi: 10.1086/527461
  * Kaastra & Bleeker (2016) Kaastra, J. S., & Bleeker, J. A. M. 2016, A&A, 587, A151, doi: 10.1051/0004-6361/201527395
  * Kennea et al. (2022) Kennea, J. A., Evans, P. A., & Negoro, H. 2022, The Astronomer’s Telegram, 15421, 1
  * King & Ritter (1998) King, A. R., & Ritter, H. 1998, MNRAS, 293, L42, doi: 10.1046/j.1365-8711.1998.01295.x
  * Kong (2022) Kong, A. K. H. 2022, The Astronomer’s Telegram, 15501, 1
  * Kuiper & Hermsen (2009) Kuiper, L., & Hermsen, W. 2009, A&A, 501, 1031, doi: 10.1051/0004-6361/200811580
  * Kuiper et al. (2020) Kuiper, L., Tsygankov, S. S., Falanga, M., et al. 2020, A&A, 641, A37, doi: 10.1051/0004-6361/202037812
  * Li et al. (2018) Li, Z., De Falco, V., Falanga, M., et al. 2018, A&A, 620, A114, doi: 10.1051/0004-6361/201833857
  * Li et al. (2022) Li, Z., Kuiper, L., Falanga, M., et al. 2022, The Astronomer’s Telegram, 15471, 1
  * Li et al. (2021) Li, Z. S., Kuiper, L., Falanga, M., et al. 2021, A&A, 649, A76, doi: 10.1051/0004-6361/202140360
  * Liu et al. (2020) Liu, C., Zhang, Y., Li, X., et al. 2020, Science China Physics, Mechanics, and Astronomy, 63, 249503, doi: 10.1007/s11433-019-1486-x
  * Mandal et al. (2023) Mandal, M., Pal, S., Chauhan, J., Lohfink, A., & Bharali, P. 2023, MNRAS, 521, 881, doi: 10.1093/mnras/stad604
  * Negoro et al. (2022) Negoro, H., Serino, M., Iwakiri, W., et al. 2022, The Astronomer’s Telegram, 15418, 1
  * Pan et al. (2022) Pan, Y. Y., Li, Z. S., Zhang, C. M., & Zhong, J. X. 2022, MNRAS, 513, 6219, doi: 10.1093/mnras/stac1365
  * Papitto et al. (2020) Papitto, A., Falanga, M., Hermsen, W., et al. 2020, New Astronomy Reviews, 91, 101544
  * Papitto et al. (2011) Papitto, A., Riggio, A., Burderi, L., et al. 2011, A&A, 528, A55, doi: 10.1051/0004-6361/201014837
  * Papitto et al. (2013) Papitto, A., Ferrigno, C., Bozzo, E., et al. 2013, Nature, 501, 517, doi: 10.1038/nature12470
  * Patruno (2010) Patruno, A. 2010, ApJ, 722, 909, doi: 10.1088/0004-637X/722/1/909
  * Patruno et al. (2009) Patruno, A., Altamirano, D., Hessels, J. W. T., et al. 2009, ApJ, 690, 1856, doi: 10.1088/0004-637X/690/2/1856
  * Patruno & Watts (2021) Patruno, A., & Watts, A. L. 2021, in ASSL, Vol. 461, Timing Neutron Stars: Pulsations, Oscillations and Explosions, ed. T. M. Belloni, M. Méndez, & C. Zhang, 143–208, doi: 10.1007/978-3-662-62110-3_4
  * Poutanen & Gierliński (2003) Poutanen, J., & Gierliński, M. 2003, MNRAS, 343, 1301, doi: 10.1046/j.1365-8711.2003.06773.x
  * Poutanen & Svensson (1996) Poutanen, J., & Svensson, R. 1996, ApJ, 470, 249, doi: 10.1086/177865
  * Powell et al. (2007) Powell, C. R., Haswell, C. A., & Falanga, M. 2007, MNRAS, 374, 466, doi: 10.1111/j.1365-2966.2006.11144.x
  * Psaltis & Chakrabarty (1999) Psaltis, D., & Chakrabarty, D. 1999, ApJ, 521, 332, doi: 10.1086/307525
  * Radhakrishnan & Srinivasan (1982) Radhakrishnan, V., & Srinivasan, G. 1982, Current Science, 51, 1096
  * Remillard et al. (2022) Remillard, R. A., Loewenstein, M., Steiner, J. F., et al. 2022, AJ, 163, 130, doi: 10.3847/1538-3881/ac4ae6
  * Sanna et al. (2020) Sanna, A., Burderi, L., Gendreau, K. C., et al. 2020, MNRAS, 495, 1641, doi: 10.1093/mnras/staa1253
  * Shapiro & Teukolsky (1983) Shapiro, S. L., & Teukolsky, S. A. 1983, Black holes, white dwarfs, and neutron stars: the physics of compact objects (New York: Wiley)
  * Tong (2015) Tong, H. 2015, Research in Astronomy and Astrophysics, 15, 517, doi: 10.1088/1674-4527/15/4/005
  * Wijnands & van der Klis (1998) Wijnands, R., & van der Klis, M. 1998, Nature, 394, 344, doi: 10.1038/28557
  * Worley et al. (2008) Worley, A., Krastev, P. G., & Li, B.-A. 2008, ApJ, 685, 390, doi: 10.1086/589823
  * Zhang et al. (2020) Zhang, S.-N., Li, T., Lu, F., et al. 2020, Science China Physics, Mechanics, and Astronomy, 63, 249502, doi: 10.1007/s11433-019-1432-6