# Unveiling Causalities in SAR ATR: A Causal Interventional Approach for
Limited Data

Chenwei Wang, Xin Chen, You Qin, Siyi Luo,
Yulin Huang, Jifang Pei, and Jianyu Yang This work was supported by the
Sichuan Science and Technology Program under Grant
2023NSFSC1970.(_Corresponding author: Yulin Huang.)_ The authors are with the
Department of Electrical Engineering, University of Electronic Science and
Technology of China, Chengdu 611731, China (e-mail<EMAIL_ADDRESS>dbw181101@163.com).

###### Abstract

Synthetic aperture radar automatic target recognition (SAR ATR) methods fall
short with limited training data. In this letter, we propose a causal
interventional ATR method (CIATR) to formulate the problem of limited SAR data
which helps us uncover the ever-elusive causalities among the key factors in
ATR, and thus pursue the desired causal effect without changing the imaging
conditions. A structural causal model (SCM) is comprised using causal
inference to help understand how imaging conditions acts as a confounder
introducing spurious correlation when SAR data is limited. This spurious
correlation among SAR images and the predicted classes can be fundamentally
tackled with the conventional backdoor adjustments. An effective implement of
backdoor adjustments is proposed by firstly using data augmentation with
spatial-frequency domain hybrid transformation to estimate the potential
effect of varying imaging conditions on SAR images. Then, a feature
discrimination approach with hybrid similarity measurement is introduced to
measure and mitigate the structural and vector angle impacts of varying
imaging conditions on the extracted features from SAR images. Thus, our CIATR
can pursue the true causality between SAR images and the corresponding classes
even with limited SAR data. Experiments and comparisons conducted on the
moving and stationary target acquisition and recognition (MSTAR) and
OpenSARship datasets have shown the effectiveness of our method with limited
SAR data.

###### Index Terms:

SAR, ATR, Causal Graph, Interventional Training

##  I Introduction

Synthetic aperture radar (SAR) is a flexible remote sensing technology, useful
in multiple civil and military contexts, offering high-resolution imaging
regardless of time or weather [1, 2, 3]. Its key application, automatic target
recognition (ATR), has evolved over the past fifty years [4, 5, 6]. Notably,
the last decade has seen considerable improvements in ATR’s performance,
largely driven by advancements in deep learning technology [7, 8, 9, 10, 11,
12, 13, 14].

While current deep learning-based SAR ATR methods exhibit encouraging
performance, they inherently rely on extensive information drawn from a large
collection of SAR images. However, the process of generating a substantial
volume of SAR images, coupled with the requirement for accurate labelling, is
both resource-intensive and time-consuming [15, 16]. As a result, there is
often a shortage of information available for supervised training, which
ultimately impacts the effectiveness of these methods. This issue underscores
a critical disconnect between the theoretical design of ATR methods and their
practical applications, with existing methods often falling short when
deployed in real-world contexts [17, 18, 19]. This problem, termed as SAR ATR
with limited training data, has recently become a focus in research [20, 21,
22, 23].

(a) Causal Graph for ATR with limited SAR data (b) Sensitiveness of SAR images
to Imaging Conditions

(c) Our CIATR

Figure 1: Causal graph for ATR with limited SAR data. (a) $I$ is SAR images,
$IC$ is the imaging conditions, $X$ is the features affected by $IC$, and $Y$
is the predicted classes. (b) Even only the azimuth angle is changing, the
scattering characteristics of SAR images varies. (c) is our CIATR as a
solution of the decounding training for ATR with limited SAR data.

The primary challenge in ATR with limited training data is the weak
performance caused by the sensitivity of SAR images to imaging conditions, as
shown in the causal graph in Fig. 1 (a). The causal graph is constructed based
on the assumption of the causalities among the SAR images $I$, the imaging
conditions $IC$, the extracted feature $X$, and the classification $Y$. When
the imaging conditions IC are changing, the inner-class SAR images have
obvious variance of scattering characteristics. Besides, it is inevitable that
the features X contains some feature IF effected by IC. As shown in Fig. 1
(b), even only a change in the azimuth angle causes inner-class SAR images to
display distinct scattering characteristics. In the process of recognition,
the ATR method aims to model $P(Y|X)$, but the imaging conditions $IC$ act as
a confounder that is the common cause of the features via $IC\to X$ and the
classification Y via $M\to IF\to Y$. As a results, the imaging conditions
introduce spurious correlation in the process of modeling $P(Y|X)$, thus
leading to the weak performance of ATR method with limited SAR data.

Therefore, in this letter, we propose a causal interventional ATR method
(CIATR) with limited SAR data that not only fundamentally analysis the role of
imaging conditions in the recognition, but also provides a principled solution
to improve the recognition performance. Specifically, our contributions are
summarized as follows.

1) Section 2.1 introduces a Structural Causal Model (SCM) that elucidates why
imaging conditions, while negligible with ample SAR data, act as a confounder
introducing spurious correlation into the ATR model when data is limited.

2) Section 2.2 then outlines specific effective implementations using backdoor
adjustment [24], which mainly consist two steps: data augmentation with
spatial-frequency domain hybrid transformation, and the feature discrimination
with hybrid similarity measurement.

3) Thanks to the causal intervention, the CIATR achieves state-of-the-art
recognition performances on MSTAR and OpenSARship data set with different
numbers of training samples. The ablation experiments have validated the
effectiveness of the CIATR.

The rest of this letter is organized below. Section II presents the causal
graph and solution for ATR with limited SAR data. Section III verifies the
effectiveness of the proposed method with experiments, and Section IV gives
our conclusion.

Figure 2: Specific implement of backdoor adjustments in our CIATR.

##  II Causal Interventional ATR Method

This section starts with the introduction of a Structural Causal Model (SCM),
detailing the cause-and-effect relationships among variables in SAR ATR. A
core solution is proposed to mitigate the false correlations arising due to
imaging conditions, thereby enhancing recognition performance even when SAR
data is scarce.

###  II-A Structural Causal Model

To systematically analyze ATR with limited SAR data, the SCM is a directed
acyclic graph that elucidates the influence of variables of interest
${\bf{I}}$, ${\bf{X}}$, $\rm{IC}$ on the ATR model’s recognition results, $Y$.
Each arrow represents a causal relationship between two nodes. The fundamental
reasoning behind SCM is detailed below.

${\bf{I}}\to Y$. The ATR modal aims to finish the precise recognition $Y$
condition on ${\bf{I}}$, $P(Y|{\bf{I}})$. There are two paths to determine $Y$
by ${\bf{I}}$: 1) ${\bf{I}}toY$ is the direct path which means ${\bf{I}}$ has
a direct effect on $Y$. 2) ${\bf{I}}\to{\bf{X}}_{v}\to Y$ is the mediation
path which means that the features ${\bf{X}}_{v}$ extracted from $\bf{I}$ play
the mediator role in the recognition process.

${\rm{IC}}\to{\bf{I}}\to Y$. The sensitivity of SAR images ${\bf{I}}$ to the
imaging conditions ${\rm{IC}}$ leads to the scattering characteristic varying
when ${\rm{IC}}$ is changing. The imaging conditions contain many aspects, for
example, the azimuth angle and the parameters of the imaging platform. Even if
one factor of imaging conditions changes, i.e., the azimuth angle, the
scattering characteristic of SAR image varies obviously. As a result, the
imaging conditions ${\rm{IC}}$ affect the scattering characteristics of the
entire SAR image ${\bf{I}}$, thereby influencing the classification $Y$.

${\rm{IC}}\to{\bf{X}}_{v}\leftarrow{\bf{I}}$. The features ${\bf{X}}_{v}$ are
denoted as the low-dimensional representations of SAR image $\bf{I}$. 1)
${\rm{IC}}\to{\bf{X}}_{v}$. The features ${\bf{X}}_{v}$ contains the features
from the targets and the background, the imaging conditions ${\rm{IC}}$ have
an obvious effect on the background, like cluster and shadow regions. Thus,
the imaging conditions ${\rm{IC}}$ and the SAR image $\bf{I}$ have the
causality for the feature ${\bf{X}}_{v}$. When modeling the recognition
$P(Y|{\bf{I}})$, the information of $\rm{IC}$ affects not only the SAR images
but also the extracted features ${\bf{X}}_{v}$.

An ideal ATR model with limited SAR data should rely on the true causality
between $\bf{I}$ and $Y$ to achieve precise recognition. However, as mentioned
above, the conventional modeling $P(Y|\bf{I})$ fails to be ideal, because the
likelihood of $Y$ given $\bf{I}$ is not only due to ${\bf{I}}\to Y$ and
${\bf{I}}\to{\bf{X}}_{v}\to Y$, but also the spurious correlation introduced
by $\rm{IC}$ via ${\bf{D}}\to{\bf{I}}$ and ${\bf{D}}\to{\bf{X}}_{v}\to Y$.

Therefore, to obtain an ideal ATR model with limited SAR data, it is necessary
to pursue the true causality between the $I$ and $Y$ without the spurious
correlation introduced by $\rm{IC}$. Fortunately, the backdoor adjustment can
be used to implement the causal intervention $P(Y|do({\bf{I}}))$ to mitigate
the spurious correlation introduced by $IC$:

$\displaystyle P(Y|do({\bf{I}}))$
$\displaystyle=\sum_{{\bf{X}}_{v}}\sum_{\rm{IC}}P(Y|{\bf{I}},{\bf{X}}_{v},{\rm{IC}})P(X_{v}|{\bf{I}},{\rm{IC}})P({\rm{IC}})$
(1)
$\displaystyle=\sum_{{\bf{X}}_{v}}\sum_{\rm{IC}}P(Y|{\bf{I}},{\bf{X}}_{v})P(X_{v}|{\bf{I}},{\rm{IC}})P({\rm{IC}})$
(2)
$\displaystyle=\sum_{\rm{IC}}P(Y|{\bf{I}},{\rm{IC}},{\bf{X}}_{v}=g(\rm{I},\rm{IC}))P({\rm{IC}})$
(3)

Due to the rule 1 of do-Calculus [25], $\mathbf{D}$ does not affect $Y$
directly, so $P(Y|\bf{I},C,D)$ in Eq. (1) can be replaced by
$P(Y|{\bf{I}},{\bf{X}}_{v})$, yielding Eq. (2). Eq. (3) is because, in our
SCM, $X_{v}$ takes a deterministic value given by function $g(I,IC)$. These
equations above means that if the imaging conditions are observable, it is
possible to employ the physical intervention by stratifying $\rm{IC}$ to
mitigate the spurious correlation introduced by $\rm{IC}$. Stratifying
$\rm{IC}$ means producing an integrated set of ${\bf{X}}_{v})$ using every
value of $IC$ for any given SAR image $i$. Therefore, in the process of
modeling $P(Y|do({\bf{I}}))$, eliminating the influence of $\rm{IC}$ and
${\bf{X}}_{v})$ can achieve Eq. (2).

In the following sections, based on the solution above, we present the
specific effectiveness implementation to improve recognition performance with
limited SAR data.

###  II-B Interventional Augmentation and Discrimination

Fig. 2 illustrates our approach. Initially, we augment data using
transformations in image and frequency domains, simulating varied imaging
conditions, a process akin to $\rm{IC}$ accumulation in Eq. (2) [25]. Next, we
apply a hybrid similarity measure for feature discrimination, calculating the
effect of different conditions on a SAR image. Using the invariant risk
minimization (IRM) concept, we provide a loss $L_{d}$ to enable ${\bf{X}}_{v}$
accumulation in Eq. (2). By minimizing $L_{d}$, our CIATR models achieve
precise recognition with limited SAR data. The details of the data
augmentation with spatial-frequency domain hybrid transformation and the
feature discrimination with hybrid similarity measurement are as follows.

Given the limited SAR training set
${\bf{D}}^{tr}=\left\\{{\bf{D}}_{1}^{tr},{\bf{D}}_{2}^{tr},..,{\bf{D}}_{C}^{tr}\right\\}$,
where ${\bf{D}}_{i}^{tr}=\left\\{{\bf{x}}_{i1},…,{\bf{x}}_{in}\right\\}$
represents the training samples of the $i$th class, and
${\bf{x}}_{ij}\in\mathbb{R}^{h\times w}$ is the $j$th sample in
${\bf{D}}_{i}^{tr}$, with $n$ being the sample number of each class.

As illustrated in Fig. 2, the spatial-frequency transformation augmentation
comprises two aspects: random frequency mask and spatial transformation.
Firstly, each image in ${\bf{D}}^{tr}$, such as ${\bf{x}}_{ij}$, undergoes a
fast Fourier transform (FFT) to derive its spectrum. We then apply a random
mask to maximize potential imaging condition estimations. The random
resolution, $rm_{re}$, denotes the smallest resolution unit in the
augmentation, implying ${\bf{x}}_{ij}$ is split into $h\times
w/({rm}_{re})^{2}$ patches. The mask ratio, ${rm}_{ra}$, indicates the number
of zeroed patches, while the location, $rm_{l}$, represents the index of these
patches. Thus, the process of random frequency mask can be presented as
${\bf{x}}^{f}_{ij}=RFM(fft({\bf{x}}_{ij}),{rm}_{re},{rm}_{ra},{rm}_{l})$,
where ${\bf{x}}^{f}_{ij}$ is the masked version of ${\bf{x}}_{ij}$,
$RFM(\cdot)$ is the operation of random frequency mask, and $fft(\cdot)$ is
the fast Fourier transform.

Then the inverse FFT is employed to convert ${\bf{x}}^{f}_{ij}$ into the
spatial domain, and the random spatial transformation with the complete
transformation set $TF=\\{{tf}_{1},…,{tf}_{Q}\\}$, where ${tf}_{i}$ is the
$i$th transformation, $Q$ is the class number of all the transformation. There
are also two random variances. The random transformation combination $q$ is
obtained by randomly sampling from $TF$, and $M_{q}$ is a random value
following a normal distribution. Thus, the process of random spatial
transformation can be presented as
${\bf{x}}^{s}_{ij}=RST(ifft({\bf{x}}^{f}_{ij}),q,M_{q})$, where
${\bf{x}}^{s}_{ij}$ is the final version of ${\bf{x}}_{ij}$ after the
augmentation, $RST(\cdot)$ is the operation of random spatial transformation,
and $ifft(\cdot)$ is the inverse FFT.

Through the process outlined above, each sample in ${\bf{D}}^{tr}$ is
augmented with an additional random version to estimate the potential impact
of various imaging conditions on SAR samples. It’s worth noting that we’ve
introduced multiple randomness factors to estimate as many imaging conditions
as possible, ensuring the comprehensiveness of stratifying $\rm{IC}$.

The limited SAR training set can be presented as
${\bf{D}}_{a}^{tr}=\\{{\bf{D}}_{1}^{tr},{\bf{D}}_{2}^{tr},..,{\bf{D}}_{C}^{tr}\\},{\bf{D}}_{i}^{tr}=\\{{\bf{x}}_{i1},{\bf{x}}^{s}_{i1},…,{\bf{x}}_{in},{\bf{x}}^{s}_{in}\\}$.
Then, a feature discrimination method employs a structural measurement to
enable the ATR model to capture effective local features. Concurrently, it
uses a vector angle measurement to enhance the discriminability of the
extracted features. The hybrid measurement aims to fulfill the summation over
${\bf{X}}_{v}$ in Eq. (2). The process of feature discrimination consists of
two parts: hybrid measurement and loss calculation, as shown in Fig. 2. For
every pair of samples in ${\bf{D}}^{tr}$, we calculate their hybrid
measurement:

$hm({\bf{x}}_{ij},{\bf{x}}_{nm})=stm({\bf{x}}_{ij},{\bf{x}}_{nm})+vam({\bf{x}}_{ij},{\bf{x}}_{nm})$
(4)

where $hm(x,y)$ is the hybrid measurement between $x$ and $y$, $stm(x,y)$ is
the structural measurement between $x$ and $y$, and $vam(x,y)$ is the vector
angle measurement. The structural similarity index measure (SSIM) is employed
as $stm(\cdot,\cdot)$, and the cosine similarity is employed as
$vam(\cdot,\cdot)$.

Then, if ${\bf{x}}_{ij}$ and ${\bf{x}}_{nm}$ belong to the same class, i.e.,
$j=m$, $hm({\bf{x}}_{ij},{\bf{x}}_{nm})$ should be as large as possible.
Conversely, if ${\bf{x}}_{ij}$ and ${\bf{x}}_{nm}$ do not belong to the same
class, i.e., $j!=m$, $hm({\bf{x}}_{ij},{\bf{x}}_{nm})$ should be as small as
possible. Thus, the discrimination loss $L_{d}$ can be calculated based on the
triplet loss. Besides, the cross-entropy loss is employed as the basic
recognition loss. The final loss is the summation of $L_{ce}$ and $L_{d}$.

Our CIATR first proposed an SCM to analyze the reason for weak performance of
ATR with limited SAR data and provides a causal solution to pursue the true
causality between the $\bf{I}$ and $Y$ without the spurious correlation
introduced by $\rm{IC}$.

##  III Experiments

In this section, we assess the effectiveness of our method using benchmark SAR
image datasets, OpenSARship and MSTAR.

###  III-A Dataset

The OpenSARship dataset, gathered from 41 diverse Sentinel-1 images,
facilitates the development of sophisticated ship detection and classification
algorithms for challenging environments. This dataset comprises 11346 slices
from 17 SAR ship types, integrated with reliable AIS information. Experiments
utilize the GRD data, featuring a $2.0{\rm m}\times 1.5{\rm m}$ resolution and
a $10{\rm m}\times 10{\rm m}$ pixel size in Sentinel-1 IW mode. Ship
dimensions span from 92m to 399m in length and 6m to 65m in width.

The MSTAR dataset, a SAR ATR performance assessment standard, was launched by
the Defense Advanced Research Project Agency and the Air Force Research
Laboratory. Acquired via Sandia National Laboratory’s STARLOS sensor, it
includes X-band SAR images with 1-ft resolution across a 0° to 360° range.

TABLE I: Image Number of Different Targets of OpenSARship Dataset Class |  Training |  Testing |  Total
---|---|---|---
Bulk Carrier |  300 |  374 |  674
Container Ship |  300 |  710 |  1010
Tanks |  300 |  253 |  553
Cargo |  300 |  557 |  857
Fishing |  300 |  121 |  421
General Cargo |  300 |  165 |  465

###  III-B Recognition Performances and Comparisons under OpenSARship and
MSTAR

In this section, the experimetns under OpenSARship and MSTAR dataset is run
and presented.

####  III-B1 Recognition Performances and Comparisons of 3 and 6 classes under
OpenSARship

The OpenSARShip dataset comprises several ship categories, representing 90% of
the international shipping market, the most common and significant ships [22].
As per [26], experiments were conducted considering different numbers of
classes: 3 and an extended set of 6. The 3-class experiment incorporates bulk
carriers, container ships, and tanks, while the extended 6-class set also
includes cargo ships, fishing vessels, and general cargo (Table I).

TABLE II: Recognition Performance (%) of 3 Classes under Different Training Data in OpenSARship Dataset Class | Training Number in Each Class
---|---
20 | 30 | 40 | 50 | 60 | 70 | 80 | 100
Bulk Carrier | 71.16 | 65.89 | 58.74 | 68.42 | 67.16 | 74.11 | 71.79 | 78.74
Container Ship | 62.52 | 70.90 | 81.63 | 81.50 | 79.53 | 78.91 | 83.23 | 78.79
Tanker | 85.88 | 86.72 | 80.79 | 81.07 | 87.29 | 83.33 | 87.01 | 88.42
Average | 70.06 | 72.87 | 74.82 | 77.62 | 77.62 | 78.48 | 80.73 | 80.85
TABLE III: Recognition Performance (%) of 6 Classes under Different Training Data in OpenSARship Dataset Class | Training Number in Each Class
---|---
20 | 30 | 40 | 50 | 60 | 70 | 80 | 100
Bulk Carrier | 64.42 | 59.16 | 48.84 | 58.95 | 65.68 | 68.42 | 61.47 | 64.00
Container Ship | 55.86 | 69.30 | 78.30 | 67.94 | 73.37 | 75.34 | 75.46 | 82.86
Tanker | 50.28 | 50.85 | 40.96 | 55.93 | 57.91 | 45.48 | 55.65 | 51.41
Cargo | 36.09 | 36.09 | 43.63 | 42.73 | 38.96 | 38.78 | 41.11 | 43.45
Fishing | 85.12 | 87.60 | 93.39 | 86.78 | 90.08 | 95.87 | 87.60 | 93.39
General Cargo | 38.79 | 50.91 | 43.64 | 42.42 | 46.67 | 53.33 | 53.94 | 44.85
Average | 52.56 | 56.95 | 57.99 | 58.07 | 61.01 | 61.10 | 61.42 | 63.91

Tables II and III illustrate our method’s superiority in 3-class and 6-class
SAR ship image recognition tasks with training samples per class between 20
and 100. The 3-class recognition rate rises from 70.06% with 20 samples to
77.62% with 50, demonstrating effective use of additional samples. In 6-class
recognition, performance ascends from 52.55% with 20 samples to 63.91% with
100.

TABLE IV: Comparison of Performances (%) of 3 classes under OpenSARShip (The number in parentheses is the number of the training samples for each method) Methods | Number range of training images in each class
---|---
1 to 50 | 51 to 100 | 101 to 338
Supervised [26] |  | 58.24 (20)
---
62.09 (40)
65.63 (80) |  | 68.75 (120)
---
70.83 (240)
CNN[21] | 62.75 (50) | 68.52 (100) | 73.68 (200)
CNN+Matrix[21] | 72.86 (50) | 75.31 (100) | 77.22 (200)
PFGFE-Net[22] | - | - | 79.84 (338)
MetaBoost[17] | - | - | 80.81 (338)
Proposed |  | 70.06 (20)
---
74.82 (40)
| 80.73 (80)
---
80.85 (100)
85.40 (200)
TABLE V: Original Image Number of Different Depressions Class | Training | Testing
---|---|---
Number | Depression | Number | Depression
BMP2-9563 | 233 | $\text{17}{}^{\circ}$ | 195 | $\text{15}{}^{\circ}$
BRDM2-E71 | 298 | 274
BTR60-7532 | 256 | 195
BTR70-c71 | 233 | 196
D7-92 | 299 | 274
2S1-b01 | 299 | 274
T62-A51 | 299 | 273
T72-132 | 232 | 196
ZIL131-E12 | 299 | 274
ZSU234-d08 | 299 |  | 274 |

In comparison with other methods (Table IV), our method excels. It yields
rates of 70.06% and 74.81% when trained with 20 and 40 samples respectively,
outperforming the Supervised method’s 58.24% with 20 samples. Furthermore, our
method surpasses Supervised, CNN, and CNN+Matrix recognition rates with 80
samples. Thus, within the 1-50 samples band, our method outmatches state-of-
the-art techniques.

In conclusion, the results under OpenSARship dataset illustrates the
resilience and effectiveness of our method with limited samples, proving its
capability to perform SAR ATR tasks under resource-constrained conditions.

####  III-B2 Recognition Performances and Comparisons under MSTAR

In this subsection, we discuss the recognition performance of our CIATR and
compare it with other algorithms under MSTAR dataset. In the K-shot setting,
both our CIATR and other methods randomly select K images from each class in
the MSTAR dataset for training in Table V.

Table VI shows the recognition results when the training samples are ranging
from 5 to 100. Examining the recognition rate relative to the growth in sample
size, we find that the rate increases significantly from 75.05% to 86.47% as
samples double from 5 to 10. The rate continues to rise, reaching 96.70% with
25 samples and 97.32% with 40 samples, demonstrating the model’s learning
capability with more data. However, rate improvement slows down beyond 40
samples, only reaching 98.68% with 80 samples, indicating a learning
saturation. Still, with 100 samples, the rate peaks at 98.89%, showing steady
improvement. This analysis confirms the robustness of our method, even with
limited samples, for SAR ATR.

Table VII is the comparison with other state-of-the-art methods for SAR ATR
with limited data. Comparatively, the recognition performance significantly
decreases when the per-class image number reduces from all data to 20 samples.
MGAN-CNN mildly improves performance under 20 and 40 samples, while
Semisupervised enhances it under 20, 40, and 80 samples, utilizing self-
consistent augmentation and training resources. Quantitatively, our CIATR
outperforms others under any number of training images per class, particularly
under limited SAR training samples.

From the recognition performances and comparisons above, the superiority and
effectiveness of our method have been validated.

TABLE VI: Recognition Performance (%) under SOC on MSTAR Class | Training Number in Each Class
---|---
5 | 10 | 25 | 30 | 40 | 60 | 80 | 100
BMP2-9563 | 47.18 | 66.15 | 80.00 | 90.77 | 87.69 | 93.33 | 94.36 | 95.90
BRDM2-E71 | 88.32 | 94.16 | 98.18 | 98.18 | 99.27 | 100.00 | 100.00 | 99.64
BTR60-7532 | 62.05 | 58.97 | 95.90 | 93.33 | 95.38 | 96.92 | 98.46 | 97.44
BTR70-c71 | 59.69 | 81.63 | 92.86 | 92.86 | 92.35 | 96.94 | 96.94 | 95.92
D7-92 | 55.11 | 92.70 | 100.00 | 99.27 | 100.00 | 99.64 | 99.64 | 100.00
2S1-b01 | 82.85 | 83.58 | 96.72 | 97.45 | 96.35 | 96.72 | 96.72 | 98.91
T62-A51 | 89.38 | 87.91 | 99.27 | 98.53 | 99.63 | 99.27 | 100.00 | 99.63
T72-132 | 69.39 | 92.35 | 98.47 | 98.98 | 97.96 | 100.00 | 98.98 | 99.49
ZIL131-E12 | 93.43 | 93.80 | 100.00 | 99.27 | 100.00 | 100.00 | 100.00 | 100.00
ZSU234-d08 | 85.40 | 100.00 | 100.00 | 100.00 | 100.00 | 100.00 | 100.00 | 100.00
Average | 75.05 | 86.47 | 96.70 | 97.24 | 97.32 | 98.47 | 98.68 | 98.89
TABLE VII: Comparison of Performance (%) under SOC of MSTAR. Algorithms | Image Number for Each Class
---|---
20 | 40 | 80 | All data
PCA+SVM [7] | 76.43 | 87.95 | 92.48 | 94.32
ADaboost [7] | 75.68 | 86.45 | 91.45 | 93.51
LC-KSVD [7] | 78.83 | 87.39 | 93.23 | 95.13
DGM [7] | 81.11 | 88.14 | 92.85 | 96.07
DNN1 [8] | 77.86 | 86.98 | 93.04 | 95.54
DNN2 [15] | 79.39 | 87.73 | 93.76 | 96.50
CNN1 [7] | 81.80 | 88.35 | 93.88 | 97.03
CNN2 [8] | 75.88 | - | - | -
CNN+matrix [8] | 82.29 | - | - | -
GAN-CNN [7] | 84.39 | 90.13 | 94.91 | 97.53
MGAN-CNN [7] | 85.23 | 90.82 | 94.91 | 97.81
Semisupervised [26] | 92.62 | 97.11 | 98.65 | -
Ours | - | 97.32 | 98.68 | -

##  IV Conclusion

In conclusion, we have proposed the causal interventional ATR method (CIATR)
to address the challenges of SAR ATR with limited training data. Our approach
leverages causal inference and backdoor adjustments to mitigate the effects of
varying imaging conditions. The structural causal model (SCM) helps us
understand the role of imaging conditions as confounders in introducing
spurious correlations in ATR models when data is limited. By using data
augmentation with spatial-frequency domain hybrid transformation and feature
discrimination with hybrid similarity measurement, we effectively estimate and
mitigate the impacts of imaging conditions on SAR image features. The CIATR
method enables us to establish causal relationships between SAR images and
their corresponding classes, even with limited data. Experimental results on
the MSTAR and OpenSARship datasets validate the effectiveness of our approach
in improving SAR ATR performance under limited training conditions.

## References

  * [1] J. C. Curlander and R. N. McDonough, _Synthetic aperture radar_. Wiley, New York, 1991, vol. 11.
  * [2] C. Wang, J. Pei, X. Liu, Y. Huang, D. Mao, Y. Zhang, and J. Yang, “Sar target image generation method using azimuth-controllable generative adversarial network,” _IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing_ , vol. 15, pp. 9381–9397, 2022.
  * [3] C. Wang, X. Liu, J. Pei, Y. Huang, Y. Zhang, and J. Yang, “Multiview attention cnn-lstm network for sar automatic target recognition,” _IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing_ , vol. 14, pp. 12 504–12 513, 2021.
  * [4] B. Bhanu, “Automatic target recognition: State of the art survey,” _IEEE transactions on aerospace and electronic systems_ , no. 4, pp. 364–379, 1986.
  * [5] C. Wang, J. Pei, Z. Wang, Y. Huang, J. Wu, H. Yang, and J. Yang, “When deep learning meets multi-task learning in sar atr: Simultaneous target recognition and segmentation,” _Remote Sensing_ , vol. 12, no. 23, p. 3863, 2020.
  * [6] C. Wang, J. Pei, M. Li, Y. Zhang, Y. Huang, and J. Yang, “Parking information perception based on automotive millimeter wave sar,” in _2019 IEEE Radar Conference (RadarConf)_. IEEE, 2019, pp. 1–6.
  * [7] C. Zheng, X. Jiang, and X. Liu, “Semi-supervised sar atr via multi-discriminator generative adversarial network,” _IEEE Sensors Journal_ , vol. 19, no. 17, pp. 7525–7533, 2019.
  * [8] D. A. Morgan, “Deep convolutional neural networks for atr from sar imagery,” in _Algorithms for Synthetic Aperture Radar Imagery XXII_ , vol. 9475\. International Society for Optics and Photonics, 2015, p. 94750F.
  * [9] J. Chen, X. Qiu, C. Ding, and Y. Wu, “Sar image classification based on spiking neural network through spike-time dependent plasticity and gradient descent,” _ISPRS Journal of Photogrammetry and Remote Sensing_ , vol. 188, pp. 109–124, 2022. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0924271622000946
  * [10] L. Li, W. Ji, Y. Wu, M. Li, Y. Qin, L. Wei, and R. Zimmermann, “Panoptic scene graph generation with semantics-prototype learning,” 2023.
  * [11] R. Liang, Y. Yang, H. Lu, and L. Li, “Efficient temporal sentence grounding in videos with multi-teacher knowledge distillation,” 2023.
  * [12] C. Wang, J. Pei, J. Yang, X. Liu, Y. Huang, and D. Mao, “Recognition in label and discrimination in feature: A hierarchically designed lightweight method for limited data in sar atr,” _IEEE Transactions on Geoscience and Remote Sensing_ , vol. 60, pp. 1–13, 2022.
  * [13] C. Wang, S. Luo, J. Pei, X. Liu, Y. Huang, Y. Zhang, and J. Yang, “An entropy-awareness meta-learning method for sar open-set atr,” _IEEE Geoscience and Remote Sensing Letters_ , 2023.
  * [14] C. Wang, J. Pei, X. Liu, Y. Huang, and J. Yang, “A deep deformable residual learning network for sar image segmentation,” in _2021 IEEE Radar Conference (RadarConf21)_. IEEE, 2021, pp. 1–5.
  * [15] J. Ding, B. Chen, H. Liu, and M. Huang, “Convolutional neural network with data augmentation for sar target recognition,” _IEEE Geoscience and remote sensing letters_ , vol. 13, no. 3, pp. 364–368, 2016.
  * [16] C. Wang, J. Pei, S. Luo, W. Huo, Y. Huang, Y. Zhang, and J. Yang, “Sar ship target recognition via multiscale feature attention and adaptive-weighed classifier,” _IEEE Geoscience and Remote Sensing Letters_ , vol. 20, pp. 1–5, 2023.
  * [17] H. Zheng, Z. Hu, J. Liu, Y. Huang, and M. Zheng, “Metaboost: A novel heterogeneous dcnns ensemble network with two-stage filtration for sar ship classification,” _IEEE Geoscience and Remote Sensing Letters_ , 2022.
  * [18] C. Wang, Y. Huang, X. Liu, J. Pei, Y. Zhang, and J. Yang, “Global in local: A convolutional transformer for sar atr fsl,” _IEEE Geoscience and Remote Sensing Letters_ , vol. 19, pp. 1–5, 2022.
  * [19] C. Wang, J. Pei, Z. Wang, Y. Huang, and J. Yang, “Multi-view cnn-lstm neural network for sar automatic target recognition,” in _IGARSS 2020-2020 IEEE International Geoscience and Remote Sensing Symposium_. IEEE, 2020, pp. 1755–1758.
  * [20] J. Geng, W. Jiang, and X. Deng, “Multi-scale deep feature learning network with bilateral filtering for sar image classification,” _ISPRS Journal of Photogrammetry and Remote Sensing_ , vol. 167, pp. 201–213, 2020. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0924271620301945
  * [21] Y. Li, X. Li, Q. Sun, and Q. Dong, “Sar image classification using cnn embeddings and metric learning,” _IEEE Geoscience and Remote Sensing Letters_ , vol. 19, pp. 1–5, 2020.
  * [22] T. Zhang and X. Zhang, “A polarization fusion network with geometric feature embedding for sar ship classification,” _Pattern Recognition_ , vol. 123, p. 108365, 2022.
  * [23] C. Wang, X. Liu, Y. Huang, S. Luo, J. Pei, J. Yang, and D. Mao, “Semi-supervised sar atr framework with transductive auxiliary segmentation,” _Remote Sensing_ , vol. 14, no. 18, p. 4547, 2022.
  * [24] M. Glymour, J. Pearl, and N. P. Jewell, _Causal inference in statistics: A primer_. John Wiley & Sons, 2016.
  * [25] J. Pearl, “The do-calculus revisited,” _arXiv preprint arXiv:1210.4852_ , 2012\.
  * [26] C. Wang, J. Shi, Y. Zhou, X. Yang, Z. Zhou, S. Wei, and X. Zhang, “Semisupervised learning-based sar atr via self-consistent augmentation,” _IEEE Transactions on Geoscience and Remote Sensing_ , vol. 59, no. 6, pp. 4862–4873, 2020.


# The formation of blue large-amplitude pulsators from white-dwarf main-
sequence star mergers

Xianfei Zhang Institute for Frontier in Astronomy and Astrophysics, Beijing
Normal University, Beijing, 102206, China  Department of Astronomy, Beijing
Normal University, Beijing, 100875, China<EMAIL_ADDRESS>C. Simon Jeffery
Armagh Observatory and Planetarium, College Hill, Armagh BT61 9DG, UK  Jie Su
Yunnan Observatories, Chinese Academy of Sciences, Kunming 650216, China  Key
Laboratory for the Structure and Evolution of Celestial Objects, Chinese
Academy of Sciences, Kunming 650216, China  International Centre of
Supernovae, Yunnan Key Laboratory, Kunming 650216, China  Shaolan Bi Institute
for Frontier in Astronomy and Astrophysics, Beijing Normal University,
Beijing, 102206, China  Department of Astronomy, Beijing Normal University,
Beijing, 100875, China

###### Abstract

Blue large-amplitude pulsators (BLAPs) are hot low-mass stars which show
large-amplitude light variations likely due to radial oscillations driven by
iron-group opacities. Period changes provide evidence of both secular
contraction and expansion amongst the class. Various formation histories have
been proposed, but none are completely satisfactory. Zhang et al. (2017)
proposed that the merger of a helium core white dwarf with a low-mass main-
sequence star (HeWD+MS) can lead to the formation of some classes of hot
subdwarf. We have analyzed these HeWD+MS merger models in more detail. Between
helium-shell ignition and full helium-core burning, the models pass through
the volume of luminosity – gravity– temperature space occupied by BLAPs.
Periods of expansion and contraction associated with helium-shell flashes can
account for the observed rates of period change. We argue that the HeWD+MS
merger model provides at least one BLAP formation channel.

stars: binary:close—stars: abundances—stars: chemically
peculiar—star:evolution—white dwarfs

##  1 Introduction

Figure 1: Location of BLAPs and evolutionary tracks of HeWD+MS merger models
in the Hertzsprung-Russell diagram. The main sequence evolutionary tracks for
1 to 10 $\rm M_{\odot}$ shown as grey lines for metallicity $Z=0.02$ (mass
fraction). The black line indicates the evolutionary track of a
$0.250$+$0.650\,\mbox{$\rm M_{\odot}$}$ HeWD+MS merger remnant. The star
indicates the first helium shell flash (He-flash). Square regions show
locations of BLAPs and different types of pulsating hot subdwarfs, including
short-period p-mode sdB stars (sdBVr: purple), long-period g-mode sdBs (sdBVs:
pink), and the hotter sdO variables (sdOV: grey).

Pulsating stars have been identified in diverse groups across the
Hertzsprung–Russell (H–R) diagram (Kurtz, 2022; Jeffery & Saio, 2016). Using
periodic light variations, asteroseismology is used to measure basic
parameters and to infer the internal structure of pulsating stars (Huber et
al., 2011; Bedding et al., 2011; Chaplin & Miglio, 2013). Due to the large
scale-surveys, more and more pulsating stars have been discovered. Many new
classes have been identified in recent years, including the blue large-
amplitude pulsators (BLAPs) discovered by Pietrukowicz et al. (2017).

BLAPs have several interesting properties:
(1) BLAPs have periods in the range 2 – 60 min. As their name implies, their
light variations show a larger amplitude than other early-type pulsators with
similar periods, i.e., from 0.05 to 0.4 mag at optical wavelengths. In
general, it is difficult for short-period stars to drive a high amplitude;
e.g., the pulsating Ap stars have periods in the range 5 – 24 min and only
have relatively lower amplitudes from 0.001 to 0.02 mag.
(2) BLAPs have higher effective temperatures (25,000-34,000K) and gravities
($\log(\rm g\,/cm\,s^{-2})$=4.2-5.7) than classical Cepheids and RR Lyrae-type
stars but have similar saw-tooth light curves.
(3) The periods of BLAPs are not constant; a key observation is the detection
of period changes with both positive and negative signs observable on
timescales of years (Pietrukowicz et al., 2017). These imply that BLAPS are
observed in both rapid expansion and contraction phases, which evolution
models must be capable of explaining.
(4) Most BLAPs have an enrichment of surface helium ($\log{n(\rm He)/n(\rm
H)}$ in the range -2.4 to -0.4), which is similar to intermediate helium-rich
hot subdwarfs.
(5) The observed BLAPs have been divided into two groups: (a) low-gravity
BLAPs (classical BLAPs) with pulsation periods in the range 20–60 mins,
$\log(\rm g\,/cm\,s^{-2})$=4.2-4.7, and amplitudes 0.2–0.4 mag; (b) high-
gravity BLAPs, with periods in the range 2–8 mins, $\log(\rm
g\,/cm\,s^{-2})$=5.3-5.7, and amplitudes 0.05-0.2 mag. The high-gravity BLAPs
reported by Kupfer et al. (2019) have spectral properties and pulsation
periods similar to p-mode hot subdwarfs (Østensen et al., 2010). Unlike most
BLAPs, TMTS-BLAP-1 is located in the period gap between low-gravity and high-
gravity BLAPs(Lin et al., 2022). BLAP OW-BLAP-1 is also found in the ”gap”
(Ramsay et al., 2022). Thus, the gap between low-gravity and high-gravity
BLAPs may not be real.
(6) They are rare, and only HD 133729 has been clearly identified as binary
(Pigulski et al., 2022). Whether any others BLAPs are binaries is not known.

Thus, a successful model of BLAPs must explain the formation channel, the
driving mechanism of pulsation, the enrichment of helium abundance, the
relation between low-gravity and high-gravity BLAPs, period changes and space
density. In the standard stellar evolution theory, obtaining a single model to
explain all of these features is difficult.

Based on stellar structure, two principal scenarios to explain BLAPs were
proposed by Pietrukowicz et al. (2017); several authors have tried to
reproduce the pulsational properties using models consistent with one or both
scenarios. (1) The shell-hydrogen-burning model: low-mass pre-white dwarfs
($\sim 0.3\,\mbox{$\rm M_{\odot}$}$) with hydrogen-burning shells pass through
the BLAP instability zone, where the opacity of iron-group elements drives
their pulsations (Wu & Li, 2018; Byrne & Jeffery, 2018; Romero et al., 2018;
Byrne & Jeffery, 2020; Byrne et al., 2021). (2) The helium-core burning model:
BLAPs have similar effective temperatures but surface gravities lower than
associated with the extended horizontal-branch. Hence, BLAPs could be evolving
toward or away from the extended horizontal-branch, either as pre- or post-
hot subdwarfs (Wu & Li, 2018; Kupfer et al., 2019; Meng et al., 2020; Xiong et
al., 2022; Lin et al., 2022). Both models can represent some features of
BLAPs, but not all.

On the Hertzsprung-Russell diagram, BLAPs are closely associated with the area
occupied by hot subdwarfs. Hot subdwarf stars can be roughly divided into
subdwarf B (sdB), subdwarf O (sdO), and subdwarf OB (sdOB) by spectrum (Heber,
2009, 2016). Most hot subdwarfs have a nearly pure hydrogen surface. Some 10%
of hot subdwarfs have a surface helium abundance $>90\%$ by number (Drilling
et al., 2013; Luo et al., 2016; Lei et al., 2019; Luo et al., 2021; Lei et
al., 2022). The formation channel of most hot subdwarfs is well explained by
the interaction of binaries (Han et al., 2002, 2003; Zhang & Jeffery, 2012).
Between the H-rich and He-rich type of stars, a small number of hot subdwarfs
have a surface helium number fraction of $10-90\%$ and are referred to as
intermediate helium-rich (iHe-rich) hot subdwarfs(Ahmad & Jeffery, 2003;
Naslim et al., 2010; Jeffery et al., 2021). It is not known whether these
represent an intermediate state of either H-rich or He-rich subdwarfs, some
other evolution channel, or a combination of several channels. For example, it
is possible that they formed from a helium white dwarf (HeWD) merge with a
main-sequence (MS) companion (Zhang et al., 2017). BLAPS show a similar
surface abundance to the intermediate helium-rich hot subdwarfs. It is
therefore helpful to investigate the properties and formation channels of iHe-
rich subdwarfs as a possible channel for the formation of BLAPs.

One such channel is the white dwarf main sequence merger. Many short-period
detached binary systems consist of a helium white dwarf (HeWD) with a main-
sequence (MS) companion (Zorotovic et al., 2011). For example, SDSS
J121010.1+334722.9 is a cool $0.4\,\mbox{$\rm M_{\odot}$}$ HeWD with a
$0.16\,\mbox{$\rm M_{\odot}$}$ M dwarf companion in a 3 h eclipsing binary.
Owing to a combination of gravitational-wave radiation, tidal interaction and
magnetic braking, the orbital period and separation of such a binary can
decrease over time and, as a consequence, the MS star may fill its Roche lobe.
If the MS star has a low mass, $M_{\rm{MS}}\leq 0.7\,\mbox{$\rm M_{\odot}$}$,
the mass transfer is expected to be dynamically unstable and lead to a merger
(Hurley et al., 2002; Shen et al., 2009). The immediate products of these
HeWD+MS mergers are expected to be red giant branch-like (RGB-like) stars
(Hurley et al., 2002). Some such remnants are expected to ignite helium with a
low envelope mass and thus become hot subdwarfs (Clausen & Wade, 2011). Zhang
et al. (2017) found that some of the mergers result in the formation of hot
subdwarfs with intermediate helium-rich surfaces.

In the Zhang et al. (2017) HeWD+MS merger model, the star takes a few tens of
Myrs following the merger before reaching the He-burning main sequence (or
zero-age extended horizontal branch). The evolutionary tracks pass through the
region of the HR diagram occupied by BLAPs, and the surface abundance of
helium, representing a mixture of hydrogen and helium, approximately matches
that of BLAPs for which measurements exist. This paper therefore examines the
products of HeWD+MS mergers in more detail in order to assess their candidacy
as BLAPs. In § 4, we introduce the post-merger evolution models. The
comparison of theory with the observation of BLAPs is shown in section § 3.
The conclusions are in § 4.

##  2 The mergers

Figure 2: A part of the post-merger evolution of a $0.250+0.650\,\mbox{$\rm
M_{\odot}$}$ HeWD+MS system during three helium shell flashes. Left panel: The
dashed line shows evolution in the $\rm log\,T_{\rm{eff}}$–$\log g$ diagram.
Solid line indicates a selected loop of evolution. Numbers indicate stages
identified in the text. Right panel: The evolution of radius (red) and helium-
shell luminosity (black) for the same model. Figure 3: Evolutionary tracks of
seven HeWD+MS merger models that pass through parameter space where BLAPs are
found. The tracks are distinguished by colour and identified by mass in the
key. Top panel: evolution in the $\rm log\,T_{\rm{eff}}$–$\log L$ (H-R)
diagram. Bottom panel: evolution in the $\rm log\,T_{\rm{eff}}$–$\log g$
diagram. Filled dots with errors represent observed BLAPs 1, 9, 11, 14 and 15
– 23 from Table 1. The triangle and circle show the stars with previously
determined luminosity, OGLE-BLAP-009 (Meng et al., 2020) and TMTS-BLAP-1 (Lin
et al., 2022), respectively.  Figure 4: Evolution of surface helium abundance
$\log\,\rm N(He)/N(H)$ during helium shell flashes in seven possible BLAP
models as a function of $\log g$. Tracks are identified by colour and key as
in Fig.3. Filled dots with errors represent observed BLAPs 1, 9, 11, 14 and 15
– 23 in Table 1.

Zhang et al. (2017) presented 29 HeWD+MS post-merger models. Some of the
evolutionary tracks of post-merger pass through the region of the HR diagram
occupied by BLAPs. Zhang et al. (2017) have not analyzed the models for
pulsation stability. An opacity bump due to iron-group element abundances
enhanced by radiation levitation is required to drive BLAP pulsations
(e.g.,Byrne & Jeffery (2020)). It is therefore worth investigating the
pulsation properties of the Zhang et al. (2017) post-merger models with full
chemical diffusion. We used the stellar evolution code mesa version 23.05.11
(Modules for Experiments in Stellar Astrophysics; Paxton et al., 2011, 2013,
2015, 2018, 2019) to calculate similar models to Zhang et al. (2017). We have
analyzed the pulsation stability of each models using GYRE (Townsend &
Teitler, 2013).

Zhang et al. (2017) describe how, after the merger, an RGB-like star forms in
which the structure is a very degenerate He core surrounded by an extended
hydrogen envelope. Following a similar evolution to normal RGB stars, hydrogen
burns in a shell, and helium ash adds to the helium core. Meanwhile mass is
lost from the surface through a stellar wind. Fresh helium continues to be
produced during RGB-like evolution to compress and heat the He core. Once the
helium core is massive enough, the helium ignites in a shell, followed by a
series of helium flashes propagating inward toward the star center. The first
helium shell-flash is the strongest. Its position near the surface of the
degenerate He core drives a strong convection zone upwards that reaches the
surface, enriching the surface in helium forced upward from beneath to create
an iHe-sdB star. The flash forces the star to expand, and then shrink as the
shell luminosity drops again, so that the evolution track performs loops in
the HR diagram (Sweigart, 1997; Brown et al., 2001; Saio & Jeffery, 2000).
After a few Myrs, the He-burning flame reaches the center and the star
commences true core He burning, i.e., on the zero-age extended horizontal
branch. While helium burns in the center, heavier elements near the surface
diffuse downwards and finally produce an almost pure hydrogen atmosphere.
Thus, the remnants evolve to become hydrogen-rich single hot subdwarfs
(H-sdB). The entire evolutionary track might be represented as: $\rm
RGB\rightarrow\rm iHe-sdB\,(BLAP)\rightarrow\rm H-sdB\rightarrow WD$

As identified by Zhang et al. (2017), the evolution of HeWD+MS mergers can be
divided into two paths corresponding to an early hot flasher or a late hot
flasher (Lanz et al., 2004; Heber, 2009, 2016). (1) late hot flasher: the
hydrogen envelope is of very low mass, so flash-driven convection can yield a
maximum helium surface abundance $Y=0.954$ (mass fraction); (2) early hot
flasher: the hydrogen shell re-ignites after the helium shell flash, the star
expands, and initiates deep opacity-driven surface convection. The combination
of flash-driven convection followed by opacity-driven convection some helium
and other newly produced elements are dredged to the surface, yielding a
maximum surface abundance $Y=0.636$.

Analyzing each model in detail, we selected seven early hot flasher models as
possible BLAPs, by having tracks which pass through a similar volume of
gravity-temperature space and $\log{n(\rm He)/n(\rm H)}\leq 0$. These models
are, $0.250$+$0.650\,\mbox{$\rm M_{\odot}$}$, $0.275$+$0.660\,\mbox{$\rm
M_{\odot}$}$, $0.300$+$0.660\,\mbox{$\rm M_{\odot}$}$,
$0.325$+$0.660\,\mbox{$\rm M_{\odot}$}$, $0.350$+$0.650\,\mbox{$\rm
M_{\odot}$}$, $0.375$+$0.650\,\mbox{$\rm M_{\odot}$}$ and
$0.400$+$0.640\,\mbox{$\rm M_{\odot}$}$, where the first and second quantity
in each pair refers to the progenitor HeWD and MS mass, respectively. The
masses for each model as they cross the BLAP region are $0.484$, $0.492$,
$0.499$, $0.511$, $0.517$, $0.527$ and $0.543\mbox{$\rm M_{\odot}$}$. Thus,
the masses of BLAPs are in the range $0.484-0.543\mbox{$\rm M_{\odot}$}$.

Fig. 1 shows a detailed post-merger evolutionary track of an early hot flasher
model, i.e., $0.250+0.650\,\mbox{$\rm M_{\odot}$}$ remnant. The flash-driven
loops span the region where BLAPs are observed and then reach the zone of sdB
stars. Fig. 2 shows an expanded section of the same evolutionary track. Part
of flash loop is marked by points 1 to 3, which represent key phases during
the loop including (1) peak of the helium flash, (2) minimum radius, and (3)
He-shell luminosity minimum. The duration of the helium shell flash is short.
It is accompanied initially by a halt in the envelope contraction (1) and a
drop in total luminosity (1)-(2), followed by envelope expansion as heat
generated in the flash is transmitted outward on a thermal timescale (2)-(3).
Once the envelope reaches radiative equilibrium after the flash, contraction
will resume at a total luminosity slightly less than before the flash.

Table 1: List of BLAPs.

No. | name | $p({\rm min})$ | $\dot{p}/p(10^{-7}\rm yr^{-1})$ | $T_{\rm eff}$ | $\log{g}$ | $\log{n(\rm He)/n(\rm H)}$ | Ref.
---|---|---|---|---|---|---|---
1 | OGLE-BLAP-001 | 28.26 | 2.90$\pm$3.70 | 30800$\pm$500 | 4.61$\pm$0.07 |  $-0.55\pm$0.05 | (1)
2 | OGLE-BLAP-002 | 23.29 |  $-19.23\pm$8.05 | - | - | - | (1)
3 | OGLE-BLAP-003 | 28.46 | 0.82$\pm$0.32 | - | - | - | (1)
4 | OGLE-BLAP-004 | 22.36 |  $-5.03\pm$1.57 | - | - | - | (1)
5 | OGLE-BLAP-005 | 27.25 | 0.63$\pm$0.26 | - | - | - | (1)
6 | OGLE-BLAP-006 | 38.02 |  $-2.85\pm$0.31 | - | - | - | (1)
7 | OGLE-BLAP-007 | 35.18 |  $-2.40\pm$0.51 | - | - | - | (1)
8 | OGLE-BLAP-008 | 34.48 | 2.11$\pm$0.27 | - | - | - | (1)
9 | OGLE-BLAP-009 | 31.94 | 1.63$\pm$0.08 | 31800$\pm$1400 | 4.40$\pm$0.18 |  $-0.41\pm$0.13 | (1)
10 | OGLE-BLAP-010 | 32.13 | 0.44$\pm$0.21 | - | - | - | (1)
11 | OGLE-BLAP-011 | 34.87 | 6.77$\pm$8.87 | 26200$\pm$2900 | 4.20$\pm$0.20 |  $-0.45\pm$0.11 | (1)
12 | OGLE-BLAP-012 | 30.90 | 0.03$\pm$0.15 | - | - | - | (1)
13 | OGLE-BLAP-013 | 39.33 | 7.65$\pm$0.67 | - | - | - | (1)
14 | OGLE-BLAP-014 | 33.62 | 4.82$\pm$0.39 | 30900$\pm$2100 | 4.42$\pm$0.26 |  $-0.54\pm$0.16 | (1)
15 | high-gravity-BLAP-1 | 3.34 | - | 34000$\pm$500 | 5.70$\pm$0.05 | -2.1$\pm$0.2 | (2)
16 | high-gravity-BLAP-2 | 6.05 | - | 31400$\pm$600 | 5.41$\pm$0.06 | -2.2$\pm$0.3 | (2)
17 | high-gravity-BLAP-3 | 7.31 | - | 31600$\pm$600 | 5.33$\pm$0.05 | -2.0$\pm$0.2 | (2)
18 | high-gravity-BLAP-4 | 7.92 | - | 31700$\pm$500 | 5.31$\pm$0.05 | -2.4$\pm$0.4 | (2)
19 | OW-BLAP-1 | 10.8 | - | 30600$\pm$2500 | 4.67$\pm$0.25 |  $-2.1\pm$0.2 | (3)
20 | OW-BLAP-2 | 23.0 | - | 27300$\pm$1500 | 4.83$\pm$0.20 |  $-0.7\pm$0.1 | (3)
21 | OW-BLAP-3 | 28.9 | - | 29900$\pm$3500 | 4.16$\pm$0.40 |  $-0.8\pm$0.3 | (3)
22 | OW-BLAP-4 | 32.0 | - | 27300$\pm$2000 | 4.20$\pm$0.20 |  $-0.8\pm$0.2 | (3)
23 | TMTS-BLAP-1 | 18.9 | 22.3$\pm$0.9 | 31780$\pm$350 | 4.90$\pm$0.06 |  $-0.66\pm$0.05 | (4)
24 | HD 133729 | 32.27 | $-11.5$ | 29000 | 4.5 | - | (5)

(1) Pietrukowicz et al. (2017), (2) Kupfer et al. (2019), (3) Ramsay et al.
(2022), (4) Lin et al. (2022), (5) Pigulski et al. (2022).

##  3 Comparison with observation

Having found that model HeWD+MS merger remnants can become BLAPs, we compare
their properties to observed examples of such stars in more detail. Table 1
shows a sample of 24 confirmed BLAPs (Pietrukowicz et al., 2017; Kupfer et
al., 2019; Lin et al., 2022; Pigulski et al., 2022; Ramsay et al., 2022).
Since our model is for merged stars, which, unless they were originally in a
triple or higher multiplicity system, should now be single, we do not
including the binary HD 133729 for comparison (Pigulski et al., 2022). Fig. 3
shows these observed BLAPs and the tracks for the seven possible BLAP models
in both the $\rm logT_{\rm{eff}}$–$\log L$ and $\rm logT_{\rm{eff}}$–$\log g$
planes. Fig. 4 compares models and observations in the $\log g$\- surface
helium abundance plane. This figure shows that the model HeWD+MS merger
remnants can explain the surface helium abundance of BLAPs, with a strong
preference for the lowest mass models: $0.250$+$0.650\,\mbox{$\rm
M_{\odot}$}$, $0.275$+$0.660\,\mbox{$\rm M_{\odot}$}$,
$0.300$+$0.660\,\mbox{$\rm M_{\odot}$}$, $0.325$+$0.660\,\mbox{$\rm
M_{\odot}$}$ and $0.350$+$0.650\,\mbox{$\rm M_{\odot}$}$. Fig. 3 and Fig. 4
show that the spaces occupied by low-gravity and high-gravity BLAPs represent
different stages of similar tracks, which indicates a possible evolutionary
relation between both types of BLAP. The TMTS-BLAP-1 and OW-BLAP-1 could also
also associated with the evolution of merger models. At least, our theoretical
models have no “gap” between low-gravity and high-gravity BLAPs.

Figure 5: Profile of opacity, the mass fraction of Fe and Ni and work function
dW/dx of stellar interior temperature. The logarithm of the opacity is shown
by solid black line, the logarithm of the mass fraction of Fe and Ni is
indicated by the dotted line, while the value of dW/dx is indicated by the
solid pink line. The blue zone indicate the approximate temperatures of the
partial ionisation opacity peaks of iron group elements.  Figure 6: The
pulsation stability of the seven HeWD+MS merger models on the $\rm
logT_{\rm{eff}}$–$\log g$ plane. The tracks are indicated by grey dotted
lines, while the unstable modes are shown by coloured circles and and
identifed in the key. The BLAPs shown in Fig. 3 are also indicated.  Figure 7:
Evolution of 7 possible BLAP models in the period-luminosiuty (top), period-
temperature (middle) and period-gravity (bottom) planes. Tracks are identified
by colour as in Fig.3. Filled dots with errors represent observed BLAPs from
Table 1.  Figure 8: The rate of period change $r\equiv\dot{P}/P$ versus
pulsation period $P$ for the seven possible BLAP models. Tracks are identified
by colour as in Fig.3. The solid line shows the same loop through a helium
shell flash for the $0.250+0.650$ $\rm M_{\odot}$ post-merger model shown in
Fig. 2. Filled dots represent BLAPs 1–14 and 23 in Table 1.

We use the oscillation code GYRE to analyze the stability of the models. The
input model for GYRE is obtained from the MESA calculation and then uses a
non-adiabatic analysis to investigate the stability of each mode. We only
calculate the stability ($\eta$) and period (P) for radial modes ($l=0$) as
the BLAPs are known to be large amplitude pulsators. Modes with $\eta>0$ are
unstable. Fig. 5 shows the logarithm of opacity (log($\kappa$)) and the work
function ($\rm dW/dx$) as a function of interior temperature for a model
around point 2 of Fig. 2. The mass fraction of iron and nickel are also
included. The location of the peak in $\rm dW/dx$ coincides with the area of
maximum opacity, which is related to the ionization of iron and nickel around
$\rm log(T/K)=5.35$. Thus, this unstable mode is driven by the
$\kappa$-mechanism related to the opacity bump due to the accumulation of
iron-group elements in the envelope. We also check the models without
radiative levitation, which modes are always stable. These results agree with
previous studies that the modes are only excited when radiative levitation is
included (Byrne & Jeffery, 2018; Romero et al., 2018; Byrne & Jeffery, 2020).
Fig. 6 shows the pulsation stability of the seven models. The locations of
unstable models are shown by colored circles. Sections of the tracks shown in
Fig. 3 which are stable against pulsation are shown as dotted lines. The
instability region satisfactorily coincides with the region where BLAPS are
observed and indicates why BLAPS might not be observed outside this region.

BLAPs are considered to be rapidly evolving stars. Both the pulsation period
and the rates of period change are important features to compare with
observation. We calculate the period $P$ and the corresponding rates of period
change $r$ for the fundamental radial mode from our models. Comparing $P$, $M$
and $R$, the classical period mean-density relation (Eddington, 1918) for
radial-mode pulsations,

$P\approx Q_{\rm
F}/\sqrt{\frac{M}{M_{\odot}}\left(\frac{R_{\odot}}{R}\right)^{3}}\,{\rm min},$
(1)

holds with the pulsation constant for the fundamental radial mode $Q_{\rm
F}\approx 47$ min.

Fig. 7 shows the observed periods compared with models in the $P-\log
T_{\rm{eff}}$ and $P-\log g$ planes. In the $P-\log T_{\rm{eff}}$ plane, the
locations of all observed BLAPs are in the region of our proposed 7 possible
merger models. Meanwhile, the same models can explain all of the BLAPs in the
$P-\log g$ plane; six out of thirteen lie within 1$\sigma$ (68% confidence
interval), twelve lie within 2$\sigma$ (95%) and all lie within 3$\sigma$
(99%) of the model predictions.

Using the same models, we calculate the expected period change $r_{i}$ from

$r_{i}\equiv\frac{\dot{P}}{P}=\frac{\Delta P}{\Delta
t}\frac{1}{P}=\frac{P_{i+1}-P_{i}}{t_{i+1}-t_{i}}\frac{1}{P_{i+1}},$ (2)

where $i$ refers to each model along a track. Fig. 8 shows the rates of period
change in the $P-r$ plane. It is a feature of post-merger models that, like
BLAPs, they show both positive and negative values of $r$ as a consequence of
the $r-T_{\rm eff}$ loops associated with helium shell flashes. The 15 BLAPs
with measured $r$ all lie within the range covered by the theoretical values
for the 7 possible BLAP models: $\sim-10^{-5}\,\rm yr^{-1}$ to $\sim
10^{-5}\,\rm yr^{-1}$. Fig. 8 highlights the loop covering one helium-shell
flash for the $0.250+0.650$ $\rm M_{\odot}$ post-merger model already shown on
the $\log T_{\rm{eff}}-\log g$ plane in Fig. 2. This demonstrates how period
decreases (negative $\dot{P}$) as the star contracts, but briefly increases
(positive $\dot{P}$) during the shell-flash itself.

From Fig. 3 to Fig. 8, most properties of BLAPs can be reasonably well
reproduced by our merger model. Furthermore, the TMTS-BLAP-1 was identified as
the “Hertzsprung Gap” of Hot Subdwarfs (Lin et al., 2022). Here, we suggest
that TMTS-BLAP-1 also formed from a merger and is currently in the pre-sdB
stage, as shown in our models.

##  4 Conclusion

From an analysis of the HeWD+MS post-merger models, we have found that some of
the early-flasher models provide a possible channel for the origin of BLAPs.
Such mergers had previously been identified as a channel to form iHe-rich hot
subdwarfs, and that is confirmed here. Some of these models transit the
parameter space occupied by BLAPs during their helium-shell flash phase prior
to full helium core burning. At this point a star represented by these models
will become an iHe-rich sdB star and eventually, the heavier elements will
sink to leave an H-rich surface as the star approaches the He-burning main
sequence as a single H-rich sdB star.

Analysis of the distribution of post-mergers evolution tracks shows that
predictions for HeWD+MS mergers are consistent with observations of almost all
recent BLAPs in terms of surface effective temperature ($T_{\rm eff}$),
surface gravity ($\rm log\,g$), surface luminosity ($\rm log\,L$), surface
helium abundance, period ($P$), and rates of period change ($r$).
Significantly, because of the cyclic expansion and contraction of helium-
shell-flashing stars, the models predict both negative and positive values for
the period change, as observed.

BLAPS that evolve from mergers are likely to be single stars. It may be that
BLAPs, especially BLAPs in binaries, can also be formed in other channels or
formed in triple systems and leave behind a merged star in a binary(Preece et
al., 2022). Recently, the BLAP HD133729 was found to be binary. We hope that
further observations will identify whether any other BLAPs have companions.

We thank the referee for the helpful suggestions and comments that improved
the manuscript. This work is supported by the grants 12073006,
12090040,12090042,12288102,12133011 and 11833006 from the National Natural
Science Foundation of China, the Joint Research Fund in Astronomy (U2031203)
under cooperative agreement between the National Natural Science Foundation of
China (NSFC) and Chinese Academy of Sciences (CAS) and the CAS ”Light of West
China” Program. We also acknowledge the science research grant from the China
Manned Space Project with No.CMS-CSST-2021-A10. Armagh Observatory and
Planetarium is supported by a grant from the Northern Ireland Department for
Communities. X.Z. thanks Jie Lin and Tao Wu for helpful conversations.

## References

  * Ahmad & Jeffery (2003) Ahmad, A., & Jeffery, C. S. 2003, A&A, 402, 335, doi: 10.1051/0004-6361:20030233
  * Bedding et al. (2011) Bedding, T. R., Mosser, B., Huber, D., et al. 2011, Nature, 471, 608, doi: 10.1038/nature09935
  * Brown et al. (2001) Brown, T. M., Sweigart, A. V., Lanz, T., Landsman, W. B., & Hubeny, I. 2001, ApJ, 562, 368, doi: 10.1086/323862
  * Byrne & Jeffery (2018) Byrne, C. M., & Jeffery, C. S. 2018, MNRAS, 481, 3810, doi: 10.1093/mnras/sty2545
  * Byrne & Jeffery (2020) —. 2020, MNRAS, 492, 232, doi: 10.1093/mnras/stz3486
  * Byrne et al. (2021) Byrne, C. M., Stanway, E. R., & Eldridge, J. J. 2021, MNRAS, 507, 621, doi: 10.1093/mnras/stab2115
  * Chaplin & Miglio (2013) Chaplin, W. J., & Miglio, A. 2013, ARA&A, 51, 353, doi: 10.1146/annurev-astro-082812-140938
  * Clausen & Wade (2011) Clausen, D., & Wade, R. A. 2011, ApJL, 733, L42, doi: 10.1088/2041-8205/733/2/L42
  * Drilling et al. (2013) Drilling, J. S., Jeffery, C. S., Heber, U., Moehler, S., & Napiwotzki, R. 2013, A&A, 551, A31, doi: 10.1051/0004-6361/201219433
  * Eddington (1918) Eddington, A. S. 1918, MNRAS, 79, 2
  * Han et al. (2003) Han, Z., Podsiadlowski, P., Maxted, P. F. L., & Marsh, T. R. 2003, MNRAS, 341, 669, doi: 10.1046/j.1365-8711.2003.06451.x
  * Han et al. (2002) Han, Z., Podsiadlowski, P., Maxted, P. F. L., Marsh, T. R., & Ivanova, N. 2002, MNRAS, 336, 449, doi: 10.1046/j.1365-8711.2002.05752.x
  * Heber (2009) Heber, U. 2009, ARA&A, 47, 211, doi: 10.1146/annurev-astro-082708-101836
  * Heber (2016) —. 2016, PASP, 128, 082001, doi: 10.1088/1538-3873/128/966/082001
  * Huber et al. (2011) Huber, D., Bedding, T. R., Stello, D., et al. 2011, ApJ, 743, 143, doi: 10.1088/0004-637X/743/2/143
  * Hurley et al. (2002) Hurley, J. R., Tout, C. A., & Pols, O. R. 2002, MNRAS, 329, 897, doi: 10.1046/j.1365-8711.2002.05038.x
  * Jeffery et al. (2021) Jeffery, C. S., Miszalski, B., & Snowdon, E. 2021, MNRAS, 501, 623, doi: 10.1093/mnras/staa3648
  * Jeffery & Saio (2016) Jeffery, C. S., & Saio, H. 2016, MNRAS, 458, 1352, doi: 10.1093/mnras/stw388
  * Kupfer et al. (2019) Kupfer, T., Bauer, E. B., Burdge, K. B., et al. 2019, ApJL, 878, L35, doi: 10.3847/2041-8213/ab263c
  * Kurtz (2022) Kurtz, D. W. 2022, ARA&A, 60, 31, doi: 10.1146/annurev-astro-052920-094232
  * Lanz et al. (2004) Lanz, T., Brown, T. M., Sweigart, A. V., Hubeny, I., & Landsman, W. B. 2004, ApJ, 602, 342, doi: 10.1086/380904
  * Lei et al. (2019) Lei, Z., Zhao, J., Németh, P., & Zhao, G. 2019, ApJ, 881, 135, doi: 10.3847/1538-4357/ab2edc
  * Lei et al. (2022) Lei, Z., He, R., Nemeth, P., et al. 2022, arXiv e-prints, arXiv:2211.12323.  https://arxiv.org/abs/2211.12323
  * Lin et al. (2022) Lin, J., Wu, C., Wang, X., et al. 2022, Nature Astronomy, doi: 10.1038/s41550-022-01783-z
  * Luo et al. (2021) Luo, Y., Németh, P., Wang, K., Wang, X., & Han, Z. 2021, ApJS, 256, 28, doi: 10.3847/1538-4365/ac11f6
  * Luo et al. (2016) Luo, Y.-P., Németh, P., Liu, C., Deng, L.-C., & Han, Z.-W. 2016, ApJ, 818, 202, doi: 10.3847/0004-637X/818/2/202
  * Meng et al. (2020) Meng, X.-C., Han, Z.-W., Podsiadlowski, P., & Li, J. 2020, ApJ, 903, 100, doi: 10.3847/1538-4357/abbb8e
  * Naslim et al. (2010) Naslim, N., Jeffery, C. S., Ahmad, A., Behara, N. T., & Şahìn, T. 2010, MNRAS, 409, 582, doi: 10.1111/j.1365-2966.2010.17324.x
  * Østensen et al. (2010) Østensen, R. H., Oreiro, R., Solheim, J. E., et al. 2010, A&A, 513, A6, doi: 10.1051/0004-6361/200913480
  * Paxton et al. (2011) Paxton, B., Bildsten, L., Dotter, A., et al. 2011, ApJS, 192, 3, doi: 10.1088/0067-0049/192/1/3
  * Paxton et al. (2013) Paxton, B., Cantiello, M., Arras, P., et al. 2013, ApJS, 208, 4, doi: 10.1088/0067-0049/208/1/4
  * Paxton et al. (2015) Paxton, B., Marchant, P., Schwab, J., et al. 2015, ApJS, 220, 15, doi: 10.1088/0067-0049/220/1/15
  * Paxton et al. (2018) Paxton, B., Schwab, J., Bauer, E. B., et al. 2018, ApJS, 234, 34, doi: 10.3847/1538-4365/aaa5a8
  * Paxton et al. (2019) Paxton, B., Smolec, R., Schwab, J., et al. 2019, ApJS, 243, 10, doi: 10.3847/1538-4365/ab2241
  * Pietrukowicz et al. (2017) Pietrukowicz, P., Dziembowski, W. A., Latour, M., et al. 2017, Nature Astronomy, 1, 0166, doi: 10.1038/s41550-017-0166
  * Pigulski et al. (2022) Pigulski, A., Kotysz, K., & Kołaczek-Szymański, P. A. 2022, A&A, 663, A62, doi: 10.1051/0004-6361/202243293
  * Preece et al. (2022) Preece, H. P., Hamers, A. S., Battich, T., & Rajamuthukumar, A. S. 2022, MNRAS, 517, 2111, doi: 10.1093/mnras/stac2798
  * Ramsay et al. (2022) Ramsay, G., Woudt, P. A., Kupfer, T., et al. 2022, MNRAS, 513, 2215, doi: 10.1093/mnras/stac1000
  * Romero et al. (2018) Romero, A. D., Córsico, A. H., Althaus, L. G., Pelisoli, I., & Kepler, S. O. 2018, MNRAS, 477, L30, doi: 10.1093/mnrasl/sly051
  * Saio & Jeffery (2000) Saio, H., & Jeffery, C. S. 2000, MNRAS, 313, 671, doi: 10.1046/j.1365-8711.2000.03221.x
  * Shen et al. (2009) Shen, K. J., Idan, I., & Bildsten, L. 2009, ApJ, 705, 693, doi: 10.1088/0004-637X/705/1/693
  * Sweigart (1997) Sweigart, A. V. 1997, ApJL, 474, L23, doi: 10.1086/310414
  * Townsend & Teitler (2013) Townsend, R. H. D., & Teitler, S. A. 2013, MNRAS, 435, 3406, doi: 10.1093/mnras/stt1533
  * Wu & Li (2018) Wu, T., & Li, Y. 2018, MNRAS, 478, 3871, doi: 10.1093/mnras/sty1347
  * Xiong et al. (2022) Xiong, H., Casagrande, L., Chen, X., et al. 2022, arXiv e-prints, arXiv:2211.01564.  https://arxiv.org/abs/2211.01564
  * Zhang et al. (2017) Zhang, X., Hall, P. D., Jeffery, C. S., & Bi, S. 2017, ApJ, 835, 242, doi: 10.3847/1538-4357/835/2/242
  * Zhang & Jeffery (2012) Zhang, X., & Jeffery, C. S. 2012, MNRAS, 419, 452, doi: 10.1111/j.1365-2966.2011.19711.x
  * Zorotovic et al. (2011) Zorotovic, M., Schreiber, M. R., & Gänsicke, B. T. 2011, A&A, 536, A42, doi: 10.1051/0004-6361/201116626


#  Homeostatic behavioural response to COVID-19 infections returns R to a set-
point of 1

Fintan Costello
School of Computer Science,
University College Dublin

Paul Watts
Department of Theoretical Physics,
National University of Ireland Maynooth

Rita Howe
School of Public Health, Physiotherapy and Sports Science,
University College Dublin

###### Abstract

One clear aspect of behaviour in the COVID-19 pandemic has been people’s focus
on, and response to, reported or observed infection numbers in their
community. We describe a simple model of infectious disease spread in a
pandemic situation where people’s behaviour is influenced by the current risk
of infection and where this behavioural response acts homeostatically to
return infection risk to a certain preferred level. This model predicts that
the reproduction rate $R$ will be centered around a median value of 1, and
that a related measure of relative change in the number of new infections will
follow the standard Cauchy distribution. Analysis of worldwide COVID-19 data
shows that the estimated reproduction rate has a median of 1, and that this
measure of relative change calculated from reported numbers of new infections
closely follows the standard Cauchy distribution at both an overall and an
individual country level.

In epidemiological models of disease spread, infection numbers at time $t$ are
a function of disease transmissibility, incubation and recovery rates (all
fixed properties of the disease), of the proportion of infectious and
susceptible individuals in the population at time $t$ (functions of the state
at time $t-1$), and of behaviour: in particular, of the average number of
contacts individuals make with others at that time, $K_{t}$. In some models
(Bertozzi et al.,, 2020) this contact number $K_{t}$ is taken as to be
constant, giving a fixed transmission rate of $\beta$; in others $K_{t}$ (or
$\beta_{t}$) is treated as a free parameter, varying with time in a way that
is not described within the epidemiological model but instead is estimated via
fitting the model to data (Ndaïrou et al.,, 2020; IHME COVID-19 forecasting
team,, 2020) or by using mobility or contact tracing datasets (Nouvellet et
al.,, 2021; Badr et al.,, 2020; Russo et al.,, 2020).

We give a simple model of how people’s behaviours (and so contact numbers)
change over time in response to their assessment of risk of infection at that
time. In this model people’s behavioural response to infection balances the
risk of infection associated with contact against the various (economic,
social and psychological) gains associated with contact. We assume that people
can estimate their risk of infection given a certain number of contacts (a
risk that depends on infection rates in the community) and that each person
has a certain constant risk or probability of infection per day, $x$, which
they are willing to accept (whose value depends on their age, health,
financial status, and so on). Each person will set their number of contacts on
a given day so that, based on their estimate of the risk per contact, their
overall risk that day is approximately $x$ (so maximising their gains from
contact without incurring unacceptable risk). We assume that actors such as
businesses or governments will behave in a similar way, balancing risk against
gain in setting policy responses to infection. In this model people will tend
to change their behaviour so that their probability of infection varies around
$x$, reducing contacts when risk is higher than $x$ but increasing contacts
when it is below $x$. Since the population risk of infection is the average of
all individual risks, the overall probability of infection will vary over time
around some constant $X/N$ (where $X$ the sum of acceptable risk levels and
$N$ the population size), and so the expected number of new infections per day
will vary around $X$. Finally, since the reproduction rate $R$ is the number
of new infections caused by an existing infected individual, with new
infections varying around a constant this model predicts that $R$ will vary
around a median value of $1$.

Two aspects of this model may be surprising. First, it goes against the common
understanding that ‘an R above one means an outbreak is growing, and below one
means that it is shrinking’ (Adam,, 2020). In this model an $R$ below $1$ does
not necessarily mean the outbreak is shrinking: instead an $R$ below $1$ leads
to an increase in contact numbers, which can cause a subsequent increase in
new infections and in $R$. Second, this model assumes that people are able to
accurately judge the probability of infection and adjust their behaviour
appropriately, contradicting the common view that ‘In making predictions and
judgments under uncertainty, people do not appear to follow the calculus of
chance or the statistical theory of prediction. Instead they rely on a limited
number of heuristics which sometimes yield reasonable judgments and sometimes
lead to severe and systematic errors’ (Tversky and Kahneman,, 1973). This
aspect of the model is motivated by our previous work suggesting that people’s
assessment of probability do in fact follow the statistical theory of
prediction, and that observed patterns of systematic error in judgement are
caused by the regressive effects of random variation or noise (Costello and
Watts,, 2014, 2016, 2018, 2019; Howe and Costello,, 2020).

Figure 1:  Median of $R$ values on each day (points) with $2.5\%$ and $97.5\%$
percentiles for country $R$ values on each day (lines). The inset shows a
histogram of $R$ values (bin size $0.1$). There are no percentile lines before
$21$ February $2020$, because only $R$ values for China are reported before
that date.

We tested this model using data from the Our World in Data COVID hub (Ritchie
et al.,, 2020) (accessed February $22$, $2022$). This dataset gives the number
of new COVID-19 infections reported each day for $225$ countries, from the
Johns Hopkins University COVID-19 Data Repository (Dong et al.,, 2020); the
reproduction rate each day for $187$ countries, estimated using a Kalman
Filter approach (Arroyo-Marioli et al.,, 2021); and the estimated stringency
of government pandemic response each day for $173$ countries, from the Oxford
COVID-19 Government Response Tracker (Hale et al.,, 2021). The median
estimated $R$ across this dataset was $1.0$ (Fig $1$) with median estimated
$R$s for each individual country being indistinguishable from $1$ in a one-
sample t-test ($t(186)=1.77,p=0.08$, $95\%$ confidence interval for the mean:
$0.95\ldots 1.003$).111Analysis code at
https://osf.io/7jx64/?view_only=a092078d4745484aab4990a6f6fe9618.

The reproduction rate $R$ on a given day $t$ is a function of $n_{t}$, the
reported number of new cases on that day. Could this $R\approx 1$ result be an
artefact of the $n_{t}$ reporting process? One problem with COVID case numbers
is the frequent reporting of $0$ new cases: just under $25\%$ of $n_{t}$
values in the dataset were $0$, with these often indicating that no reporting
took place that day: a number of countries had reliable patterns of $n_{t}=0$
on weekend days only. These reporting gaps are visible as a spike in $R$
values at $0$ in the Fig $1$ histogram of $R$ values. To eliminate these
reporting gaps we reran our analysis on a cleaned dataset including only days
with $n_{t}>0$. The median estimated $R$ for $n_{t}>0$ was $1.02$, with median
estimated $R$s for each individual country being indistinguishable from $1$ in
a one-sample t-test ($t(186)=1.1,p=027$, $95\%$ confidence interval:
$0.99\ldots 1.03$). All subsequent analyses use this cleaned dataset.

Perhaps this result could be caused by a relationship between $R$ and the
number of tests being carried out? If testing increases when $R$ is high,
tests would include more cases likely to be negative and would reduce the
apparent value of $R$; similarly, low test numbers could increase the apparent
value of $R$. We calculated the correlation between $R$ and number of daily
tests carried out in the cleaned dataset. There was no significant difference
between median $R$ values for countries where this correlation was positive
and those where it was negative ($t(126.96)=0.36,p-=0.72$).

Perhaps this result is a consequence of government interventions alone, rather
than behavioural responses to risk? To check this we compared median $R$
values for countries where the average government stringency level was above
the overall mean stringency, and those for which it was below. There was no
significant difference between median $R$ values for these groups
($t(157.43)=1.74,p=0.08$); the average median $R$ was slightly higher in the
high-stringency group ($1.02$) than the low-stringency group ($1.0$).

These results suggest that the number of new infections at time $t$, $i_{t}$,
varies around some constant $X$ as a consequence of people’s behavioural
response to infection risk (and so $R$ is distributed around $1$). What can we
say about the distribution of these values $i_{t}$? This behavioural response
has a natural lag, $L$, which represents the time between an infection
occurring (at time $t-L$, say) and that infection being observed by others and
causing a behavioural response (at time $t$). This lag falls somewhere between
the incubation and recovery period for the infection (an infection becoming
observable only after incubation, and not being observable after recovery),
and means that the observed rate of new infections at time $t$ is equal to the
actual rate of new infections at time $t-L$. If $i_{t-L}>X$ then the overall
behavioural response at time $t$ will reduce contact numbers, pushing $i_{t}$
downwards, while if $i_{t-L}<X$ then the overall behavioural response at time
$t$ will increase contact numbers, pushing $i_{t}$ upwards, and so the
difference $i_{t}-i_{t-L}$ varies around $0$. Since this overall behavioural
response is the sum of all individual responses in the population, from the
Central Limit theorem this difference $i_{t}-i_{t-L}$ will follow a Normal
distribution $i_{t}-i_{t-L}\sim\mathcal{N}(0,\sigma_{t}^{2})$ with some
variance $\sigma^{2}_{t}$ (which may change over time). The difference
$i_{t-2L}-i_{t-L}$ will follow the same distribution (albeit with variance
$\sigma_{t-L}^{2}$).

Defining a measure of relative change in new infection numbers from time
$t-2L$ to time $t$,

$D_{L}(t)=\frac{(i_{t}-i_{t+L})-(i_{t-2L}-i_{t-L})}{(i_{t}-i_{t-L})+(i_{t-2L}-i_{t-L})}=\frac{i_{t}-i_{t-2L}}{i_{t}+i_{t-2L}-2i_{t-L}}$

we see that $D_{L}$ is the ratio of two standard Normal variables (sums of
common standard deviations cancelling) and so this measure $D_{L}$ will follow
the standard Cauchy distribution $C$ (with location $0$ and scale $1$) for $L$
between the incubation and recovery times. Assuming that changes in $n_{t}$
are proportional to changes in $i_{t}$, we predict that $D_{L}$ values
calculated from reported number of positive tests $n_{t}$ will also follow the
standard Cauchy distribution $C$.

To test this prediction we compare $D_{L}$ values calculated for our cleaned
dataset against the theoretical distribution $C$. For each country, at each
day $t$ we calculated $D_{L}(t)$ for various values of $L$. For some days
$D_{L}(t)$ could not be calculated (because one of the component infection
numbers was missing), or involved division by $0$; these values of $D_{L}(t)$
were dropped from analysis. Figure $2$ (inset) shows a probability-probability
plot comparing the cumulative probability of $D_{L}$ for $L=6$ against that of
$C$. Correlation of cumulative probabilities is a measure of goodness of fit
between observed and theoretical values; here the correlation was high
($r=0.9997$). Since probability-probability plots overweight extreme values,
we also analysed the relationship between $C$ and $D_{L}$ for values near the
midpoint of the range, by selecting the subset of $D_{L}$ values between $-15$
and $15$ (over $95\%$ of the total sample). Figure $2$ (main) shows a
histogram of these values. The correlation between $D_{L}$ and $C$ values for
this central-region histogram was $r=0.993$. As an additional check we
calculated location and scale estimates by taking the median of $D_{L}$ and
the median of the absolute value of $D_{L}$; these values were $0.02$ and
$1.02$ respectively, confirming the fit to the standard Cauchy distribution.

Figure 2:  Histogram of $D_{6}$ calculated from the cleaned dataset in the
central $-15\ldots 15$ range (bin size $0.5$) with standard Cauchy
distribution $C$ (dashed line, $C$ distribution scaled by bin size and total
histogram frequency for comparison). The inset shows a probability-probability
plot comparing theoretical and observed cumulative probabilities across the
entire range: the solid line in that plot is actually made up of $88,832$
points, one for each $D_{L}$ value calculated in the dataset: the dashed
diagonal line (mostly hidden by these points) is the line of identity between
theoretical and observed cumulative probabilities.

We carried out the same analysis for individual country data for $L=6$.
Correlations between cumulative probabilities of $C$ and $D_{L}$ for
individual countries were greater than $r=0.98$ for all countries; for binned
$D_{L}$ and $C$ values in the $-15\ldots 15$ region the mean correlation was
$r=0.97$ (with a $2.5\%\ldots 97.5\%$ percentile range of $0.83$ to $0.99$).
The $2.5\%\ldots 97.5\%$ percentile range for $D_{L}$ medians (location
parameter estimates) across countries was $-0.06\ldots 0.22$ with an overall
median of $0.03$, for absolute medians (scale parameter estimates) the same
range was $0.71\ldots 1.38$ with an overall median of $1$, agreeing with the
predicted values of $0$ and $1$ for location and scale. The standard Cauchy
distribution gives a good fit to values of $D_{L}$ for individual country
data.

We predicted that this agreement should hold for $L$ approximately between $5$
and $14$ (estimates for the incubation and recovery period for COVID-19). To
test this we carried out the above analysis for all values of $L$ from $2$ to
$100$. Probability-probability correlations, slopes and intercepts did not
change noticeably because changes in $L$ primarily affect $D_{L}$ around the
median. Figures $3$ plots the correlation between binned values of $D_{L}$ and
$C$ values in the $-15\ldots 15$ median region, for each value of $L$.
Correlation was highest for values of $L$ in the $5$ to $14$ region, as
predicted, supporting the behavioural response model.

Figure 3:  Each point shows the correlation $r$ between values of $D_{L}$ in
the $-15\ldots 15$ central range and the standard Cauchy distribution density
for those bins, for values of $L$ in from $2$ to $100$.

There are a number of clear limitations to these results. First, our model of
homeostasis due to behavioural response assumes that the susceptible
population is aware of and responding to the risk of infection, and so applies
to epidemic or pandemic situations only: we do not expect this homeostatic
effect to hold in narrower outbreak situations. Second, our model depends on
the assumption that new infection numbers at time $t$ are a reflection of the
probability of infection at that time. This assumption holds for infections
with short incubation and recovery periods; for infections where these periods
are longer, this assumption doesn’t hold. Third: our model assumes that people
are free to limit their number of contacts to match their acceptable level of
risk. For some demographics this is not the case: people in poverty, for
example, may be economically unable to limit their contacts in this way, and
so will have an estimated risk of infection systematically above their
acceptable risk level. Assuming that people’s acceptable risk levels are well-
calibrated, this predicts increased infections in such demographics (Patel et
al.,, 2020; Little et al.,, 2021). Fourth: we assume that reported infection
numbers are proportional to actual infection rates. If reported infection
numbers do not follow actual infection numbers, we do not expect these results
to hold.

A final caveat concerning the interpretation of these results. At first glance
our results may suggest that government responses to infection have no value
or no effect. This is not the case: government restrictions on contact clearly
act to reduce the risk of infection. Instead, our model suggests that
government restrictions act to reduce the overall level of acceptable risk
$X$, so reducing the number of infections to a lower (but still approximately
constant) value. Relaxation of those restrictions then produces an increase in
the overall level of acceptable risk, $X$, and so causes infection rates to
rise, but only to that new level.

## References

  * Adam, (2020) Adam, D. (2020).  A guide to r–the pandemic’s misunderstood metric.  Nature, 583(7816):346–349.
  * Arroyo-Marioli et al., (2021) Arroyo-Marioli, F., Bullano, F., Kucinskas, S., and Rondón-Moreno, C. (2021).  Tracking r of covid-19: A new real-time estimation using the kalman filter.  PloS one, 16(1):e0244474.
  * Badr et al., (2020) Badr, H. S., Du, H., Marshall, M., Dong, E., Squire, M. M., and Gardner, L. M. (2020).  Association between mobility patterns and covid-19 transmission in the usa: a mathematical modelling study.  The Lancet Infectious Diseases, 20(11):1247–1254.
  * Bertozzi et al., (2020) Bertozzi, A. L., Franco, E., Mohler, G., Short, M. B., and Sledge, D. (2020).  The challenges of modeling and forecasting the spread of covid-19.  Proceedings of the National Academy of Sciences, 117(29):16732–16738.
  * Costello and Watts, (2014) Costello, F. and Watts, P. (2014).  Surprisingly rational: Probability theory plus noise explains biases in judgment.  Psychological Review, 121(3):463–480.
  * Costello and Watts, (2016) Costello, F. and Watts, P. (2016).  People’s conditional probability judgments follow probability theory (plus noise).  Cognitive Psychology, 89:106–133.
  * Costello and Watts, (2018) Costello, F. and Watts, P. (2018).  Invariants in probabilistic reasoning.  Cognitive Psychology, 100:1–16.
  * Costello and Watts, (2019) Costello, F. and Watts, P. (2019).  The rationality of illusory correlation.  Psychological review, 126(3):437.
  * Dong et al., (2020) Dong, E., Du, H., and Gardner, L. (2020).  An interactive web-based dashboard to track covid-19 in real time.  The Lancet infectious diseases, 20(5):533–534.
  * Hale et al., (2021) Hale, T., Angrist, N., Goldszmidt, R., Kira, B., Petherick, A., Phillips, T., Webster, S., Cameron-Blake, E., Hallas, L., Majumdar, S., et al. (2021).  A global panel database of pandemic policies (oxford covid-19 government response tracker).  Nature Human Behaviour, 5(4):529–538.
  * Howe and Costello, (2020) Howe, R. and Costello, F. (2020).  Random variation and systematic biases in probability estimation.  Cognitive Psychology, 123:101306.
  * IHME COVID-19 forecasting team, (2020) IHME COVID-19 forecasting team (2020).  Modeling covid-19 scenarios for the united states.  Nature medicine.
  * Little et al., (2021) Little, C., Alsen, M., Barlow, J., Naymagon, L., Tremblay, D., Genden, E., Trosman, S., Iavicoli, L., and van Gerwen, M. (2021).  The impact of socioeconomic status on the clinical outcomes of covid-19; a retrospective cohort study.  Journal of community health, 46(4):794–802.
  * Ndaïrou et al., (2020) Ndaïrou, F., Area, I., Nieto, J. J., and Torres, D. F. (2020).  Mathematical modeling of covid-19 transmission dynamics with a case study of wuhan.  Chaos, Solitons & Fractals, 135:109846.
  * Nouvellet et al., (2021) Nouvellet, P., Bhatia, S., Cori, A., Ainslie, K. E., Baguelin, M., Bhatt, S., Boonyasiri, A., Brazeau, N. F., Cattarino, L., Cooper, L. V., et al. (2021).  Reduction in mobility and covid-19 transmission.  Nature communications, 12(1):1–9.
  * Patel et al., (2020) Patel, J., Nielsen, F., Badiani, A., Assi, S., Unadkat, V., Patel, B., Ravindrane, R., and Wardle, H. (2020).  Poverty, inequality and covid-19: the forgotten vulnerable.  Public health, 183:110.
  * Ritchie et al., (2020) Ritchie, H., Mathieu, E., Rodes-Guirao, L., Appel, C., Giattino, C., Ortiz-Ospina, E., Hasell, J., Macdonald, B., Beltekian, D., and Roser, M. (2020).  Coronavirus pandemic (covid-19).  Our World in Data.  https://ourworldindata.org/coronavirus.
  * Russo et al., (2020) Russo, L., Anastassopoulou, C., Tsakris, A., Bifulco, G. N., Campana, E. F., Toraldo, G., and Siettos, C. (2020).  Tracing day-zero and forecasting the covid-19 outbreak in lombardy, italy: A compartmental modelling and numerical optimization approach.  Plos one, 15(10):240–649.
  * Tversky and Kahneman, (1973) Tversky, A. and Kahneman, D. (1973).  Availability: A heuristic for judging frequency and probability.  Cognitive Psychology, 5:207–232.