# Image Classification with Rotation-Invariant Variational Quantum Circuits

Paul San Sebastian Sein  Ikerlan Technology Research Centre, Basque Research
and Technology Alliance (BRTA), Paseo J.M. Arizmediarrieta 2, E-20500,
Arrasate-Mondragón, Spain  Mikel Cañizo  Ikerlan Technology Research Centre,
Basque Research and Technology Alliance (BRTA), Paseo J.M. Arizmediarrieta 2,
E-20500, Arrasate-Mondragón, Spain  Román Orús  Donostia International Physics
Center, Paseo Manuel de Lardizabal 4, E-20018 San Sebastián, Spain  Multiverse
Computing, Paseo de Miramón 170, E-20014 San Sebastián, Spain  Ikerbasque
Foundation for Science, Maria Diaz de Haro 3, E-48013 Bilbao, Spain

###### Abstract

Variational quantum algorithms are gaining attention as an early application
of Noisy Intermediate-Scale Quantum (NISQ) devices. One of the main problems
of variational methods lies in the phenomenon of Barren Plateaus, present in
the optimization of variational parameters. Adding geometric inductive bias to
the quantum models has been proposed as a potential solution to mitigate this
problem, leading to a new field called Geometric Quantum Machine Learning. In
this work, an equivariant architecture for variational quantum classifiers is
introduced to create a label-invariant model for image classification with
$C_{4}$ rotational label symmetry. The equivariant circuit is benchmarked
against two different architectures, and it is experimentally observed that
the geometric approach boosts the model’s performance. Finally, a classical
equivariant convolution operation is proposed to extend the quantum model for
the processing of larger images, employing the resources available in NISQ
devices.

††preprint: APS/123-QED

##  I Introduction

Together with the recent surge of quantum technologies, both in theoretical
frameworks and practical applications, quantum computing has been placed in
the spotlight. Many quantum algorithms have shown significant speed-ups with
respect to classical algorithms, and a wide variety of applications are being
explored, from simulation of many-body systems to combinatorial optimization,
cryptanalysis, database search, machine learning, finance, and beyond [1, 2,
3, 4, 5, 6, 7, 8]. However, many of these algorithms require hardware
resources that stand far from the capabilities of the currently available
noisy intermediate-scale quantum (NISQ) devices. For this reason, and in
parallel to foundational algorithmic research focused on a quantum advantage
in fault-tolerant quantum devices, a significant effort is put into the search
for mid-term applicability and utility of algorithms in the NISQ era. Although
its computational advantage is yet to be proven, Quantum Machine Learning
(QML) seems one of the most promising candidates for an early application of
quantum devices [9, 10, 11], since in many algorithms the goal is to
approximate a distribution and the result often does not rely on an exact
value of the computation. This fact makes QML algorithms resilient to noise,
and some signs hint that the noise could even improve the optimization in
specific cases [12, 13].

QML encompasses a broad meaning, but the scope of this work will focus on the
employment of quantum circuits in learning problems that involve classical
data. Variational quantum algorithms [14, 15, 16, 17] hold a renowned place in
this paradigm, where quantum circuits are used to create expressive high-
dimensional feature maps that have been proved to be asymptotically universal
approximators [18, 19]. However, these variational algorithms suffer
trainability issues called Barren Plateaus (BP), where the variance of the
loss function tends to vanish exponentially with the size of the quantum
circuit [20, 21]. Extensive work has been done characterizing this phenomenon,
identifying its sources [22, 23, 24, 25, 26] and looking for methods to
mitigate or avoid it [27, 28, 29].

Geometric Quantum Machine Learning (GQML) [30, 31, 32, 33], has shown a great
potential to alleviate Barren Plateaus and bring significant improvements in
trainability and generalization of Variational Quantum Classifiers [34].
Classical Geometric Deep Learning [35] founded an algebraic framework to
analyze learning problems and proposed a geometric approach to add inductive
bias [36] to learning models, i.e., a blueprint to design algorithms based on
the prior knowledge about the geometric properties of the data set. As its
quantum counterpart, GQML has established a potential research line where many
authors implemented the geometric approach to a variety of practical problems,
showing promising results. For example, Ref.[37] studied the behavior of
equivariant quantum neural networks (EQNN) under the presence of noise and
suggested a strategy to enhance the protection of symmetry; Ref.[38]
introduced an EQNN for a classification task with $\mathbb{Z}_{2}$ permutation
symmetry; in Ref.[39] Quantum Fourier Transform was employed to build a
rotation equivariant QNN for scanning tunneling microscope images; Ref.[40]
built a $SU(2)$\- equivariant quantum circuit to learn on spin networks;
Ref.[41] proposed a data-encoding strategy for Point Clouds, which is
invariant under point permutations; and Ref.[42] designed a permutation-
equivariant quantum circuit to solve the Traveling Salesman Problem and train
it using Reinforcement Learning. Related to this work, Ref.[43] presented an
equivariant quantum classifier, in this case, for images with reflection
label-symmetry. Additionally, Ref.[44] covered the case of images with $p4m$
label-symmetry.

In this context, the aim of this paper is to introduce an architecture for a
variational quantum circuit for image classification, equivariant to $C_{4}$
image-rotation label-symmetry. In difference to the work of Ref.[43] and
Ref.[44], which both of them used Ampitude Encoding, this work employs a
different encoding strategy, namely Angle Encoding as discussed in III.2. This
encoding allows the application of re-uploading layers, boosting the
approximation ability of the model. Moreover, reflection symmetry is not
included in our scope, making the proposed model distinguish chiral images.
This makes the model fit, for example, for medical vision problems like brain
tumor image detection or fingerprint classification, where it is vital to
discern between the left and right sides of the body. Concretely, the
contributions of this work are as follows:

  * •

The design of an $C_{4}$ invariant Variational Quantum Classifier. A detailed
explanation of the followed methodology, as well as a benchmarking and a
comparison with the other two architectures.

  * •

A $C_{4}$ equivariant convolution operation to extend the label-invariant
quantum model for large-scale multichannel images.

The paper is organized as follows: in Sec.II, a formal explanation of the
classification problem is given, the general scheme of Variational Quantum
Classifiers (VQC) is presented, and the intuition behind GQML is summarized.
In Sec.III, the concrete application of the methodology is explained, and the
benchmarking results are shown. First, in Sec.III.1, the characteristics of
the synthetic data set generated to design the circuit are presented. Then, in
Secs.III.2 and III.3, the chosen encoding strategy and the design of the
equivariant VQC are detailed. In Sec.III.4 the comparison between results of
different models run in simulators is done, and in Sec.III.5, a method to
extend the model to process larger images is proposed. Finally, in Sec.III,
some conclusions and suggestions about future work are given.

##  II Preliminaries

###  II.1 Classification problem

Machine Learning models are meant to emulate a decision-making process that
follow intelligent beings. In this process, a response $y\in\mathcal{Y}$ is
computed as a function of a piece of input information $x\in\mathcal{X}$.
Therefore, there exists a function, often called hypothesis, that computes
$h:\;\mathcal{X}\longrightarrow\mathcal{Y}$ mapping. A model is an algorithm
designed to give a response $\hat{f}\in\mathcal{Y}$, which does not have to
coincide necessarily with $y$ in an initial stage, for an input data-point
$x$. A key characteristic of these algorithms is that they contain intern
parameters $\theta$ that can be tuned to modify their output for a given data
point. Hence, the mapping $\hat{h}:\;\mathcal{X}\longrightarrow\mathcal{Y}$
that it computes can be modified. The main task in Machine Learning is to find
the optimal configuration of $\theta$ parameters to make the model’s output
the most similar possible to the output of the target process for all inputs;
or, in other words, to make $\hat{h}$ approximate $h$. This is what training
the model means, or what it is referred to when stating that a machine can
learn.

In supervised learning, a data set
$\mathcal{D}=\\{(x^{(1)},y^{(1)}),\dots,(x^{(N)},y^{(N)})\\}$ is provided with
$N$ sample data-points and their correspondent responses. A loss function
$\mathcal{L}$ is defined, which quantifies the dissimilarity between the
hypothesis function and the approximation of the model based on the available
data set. The training of the model can be translated into an optimization
problem where the optimal parameters will be the ones that minimize the loss
function, following the so-called empirical risk minimization strategy[45]:

$\theta^{*}=arg\,\min_{\theta}\sum_{i}^{N}\mathcal{L}(y^{(i)},\hat{y}(x^{(i)},\theta)).$
(1)

In classification problems, the response is often discretized into categorical
(non-numerical) labels. The most straightforward case is the binary
classification problem, where the continuous output $f(x,\theta)$ of an
algorithm can be discretized, establishing a decision threshold, for example,
$\hat{y}(x,\theta)=sign(f(x,\theta))$.

###  II.2 Variational Quantum Classifier

Variational Quantum Classifiers (VQC) are hybrid algorithms, which make use of
quantum machines to give a prediction out from classical data and employ
classical computers to optimize the tuneable parameters. On the one hand, the
quantum circuit is composed of an encoding unitary $U(x)$ that encodes
classical data features into a quantum state. On the other hand, a set of
parameterized gates is applied, preceded by multi-qubit gates that generate
entangling, altogether forming a variational unitary $W(\theta)$ [11, 15]. The
encoding and variational unitaries will evolve the initial state of the system
into another quantum state, function of the classical input data and the
variational parameters:

$|\psi(x,\theta)\rangle=W(\theta)U(x)|0\rangle^{\otimes n_{q}},$ (2)

where $n_{q}$ is the number of qubits in the quantum circuit. The predicted
output of a data-point is defined as the expected value of an observable $O$
for the quantum state $|\psi(x,\theta)\rangle$:

$\displaystyle f(x,\theta)$ $\displaystyle=$
$\displaystyle\langle\psi(x,\theta)|O|\psi(x,\theta)\rangle$ (3)
$\displaystyle=$ $\displaystyle\langle
0|U^{\dagger}(x)W^{\dagger}(\theta)OW(\theta)U(x)|0\rangle.$

The outputs of the quantum circuit are then evaluated in the loss function,
where a classical optimizer will determine the update of variational
parameters, making use of techniques such as parameter shift [46] to compute
the gradient of the loss function respect to the optimizable parameters.

Repeating the sequence of encoding and variational unitaries, widely known as
data re-uploading [18] (see Fig. 1), has shown a significant improvement in
the model’s approximation capability. Furthermore, [19] claim that a data re-
uploading model can be considered as a Fourier series, where the accessible
frequency spectrum is directly related to the number of re-uploading layers.
Therefore, data re-uploading VQCs can be considered asymptotically universal
approximators.

Figure 1:  A general scheme of a data re-uploading VQC of $4$ qubits and $l$
layers, where a layer is defined as a sequence of data encoding and
variational unitaries.

###  II.3 Geometric Quantum Machine Learning

It is said that a data set contains a symmetry if a certain property of the
data set is preserved under the action of an operation, namely a symmetry
operation. The set of all such operations $g$ along with a composition
operation between them, constitute a group $G$ [47]. It is important to keep
in mind that this symmetry group $G$ is defined in general and abstract terms,
and it can act over different mathematical structures through its different
representations. To illustrate this, it can be thought of two representations
of the same data-point linked by a data encoding unitary $U(x)$. This unitary
maps a data-point from its original space $\mathcal{X}$ to the Hilbert space
$\mathcal{H}$ of our quantum circuit. The symmetry operation $g\in G$ then
will have two representations acting on each one of the spaces the data-point
is represented:

$\begin{CD}\mathcal{X}@>{V_{g}}>{}>\mathcal{X}\\\
@V{}V{U(x)}V@V{}V{U(x^{\prime})}V\\\
\mathcal{H}@>{U_{g}}>{}>\mathcal{H}\end{CD}$ (4)

where $V_{g}$ is the representation of $g\in G$ acting on the original data-
space and $U_{g}$ its representation acting on the Hilbert’s space. Quantum
gates are represented as unitary square matrices from the group $U(2^{n})$,
where $n$ stands for the number of qubits they act on. In consequence, it is
possible to represent $U_{g}$ as a $2^{n}\times 2^{n}$ unitary matrix, but it
is remarkable how its ultimate representation will be determined by the choice
of the data encoding method $U(x)$.

The purpose of this framework is to provide some geometric tools to design
label-invariant VQCs respect to the action of $G$:

$\displaystyle y(x)=y(V_{g}[x])$ $\displaystyle\Longrightarrow$ $\displaystyle
f(x,\theta)=f(V_{g}[x],\theta)$ (5) $\displaystyle\forall
x\in\mathcal{X}\;\textrm{and}\;\forall g\in G,$

being $f$ the output response of the model. To achieve this goal, the first
step is to ensure that the chosen encoding unitary is equivariant respect to
the symmetry group [32]:

$\displaystyle U(V_{g}[x])=U_{g}U(x)U_{g}^{\dagger}$ $\displaystyle\forall
x\in\mathcal{X}\;\textrm{and}\;\forall g\in G.$ (6)

If the condition of eq. II.3 is satisfied, an equivariant re-uploading model
can be built, symmetrizing all the gates used in the variational part. This is
done using the _Twirling formula_ [32, 30]:

$\mathcal{T}[A]=\frac{1}{|G|}\sum_{g\in G}U_{g}AU_{g}^{\dagger},$ (7)

which, applied to a unitary $A$, returns another unitary that commutes with
all symmetry operations in the symmetry group
$[\mathcal{T}(A),U_{g}]=0\;\;\forall g\in G$. Accordingly, if the variational
unitary $W$ is built out from symmterized gates, it will satisfy
$[W(\theta),U_{g}]=0$, and all re-uploading layers
$U_{l}(x,\theta_{l})=U(x)W_{l}(\theta_{l})$ will be equivariant respect the
symmetry group.

Finally, if an invariant initial-state $|0\rangle=U_{g}|0\rangle$ and an
invariant observable $O=U_{g}OU_{g}^{\dagger}$ are chosen, both for $\forall
g\in G$, a label-invariant VQC will be obtained:

$\displaystyle f(V_{g}[x],\theta)$ $\displaystyle=$ $\displaystyle\langle
0|U^{\dagger}(V[x])W^{\dagger}(\theta)OW(\theta)U(V[x])|0\rangle$ (8)
$\displaystyle=$ $\displaystyle\langle
0|U_{g}U^{\dagger}(x)U_{g}^{\dagger}W^{\dagger}(\theta)OW(\theta)U_{g}U(x)U_{g}^{\dagger}|0\rangle$
$\displaystyle=$
$\displaystyle\langle\psi(x,\theta)|U_{g}^{\dagger}OU_{g}|\psi(x,\theta)\rangle$
$\displaystyle=$ $\displaystyle\langle\psi(x,\theta)|O|\psi(x,\theta)\rangle$
$\displaystyle=$ $\displaystyle f(x,\theta).$

Recall that as $U_{g}$ is unitary, $U_{g}^{\dagger}=U_{g}^{-1}$ is also a
symmetry operation and $U_{g}^{-1}=U_{g^{-1}}$ . This property has been used
to proof the label-invariance in eq. 8.

##  III Methodology

The objective of this work is to apply the tools explained in Sec. II.3 to an
image classification problem with 4-fold rotational label-symmetry and observe
experimentally the existence of an improvement in the performance of symmetry-
preserving models. As it will be discussed later in the article, the number of
pixels encoded in the quantum circuit is a limiting factor for the simulation
of quantum circuits in classical computers. For this reason, the size of the
problem has been constrained to a reduced dimension, but the obtained results
can be generalized to larger image resolutions without loss of generality.

###  III.1 Synthetic toy dataset

The employed data set will be formed by square images of $n\times n$ pixels,
where each pixel will take values laying between 0 and 255 in grayscale.
Hence, each image can be represented by a matrix $x\in GL_{n}(\mathbb{R})$.
For a first stage, images with resolution $n=4$ of T and L tetrominoes have
been generated (see fig.2). Each tetromino can be placed in 6 different
positions and 4 orientations inside an image of $4\times 4$ pixels, but it is
possible to augment the size of the data set arbitrarily by adding a layer of
noise to each image.

(a) L tetromino

(b) T tetromino

Figure 2: Two samples of the data set.

Tetrominoes contain a 4-fold rotation label-symmetry, but notice how
reflection symmetry is not fulfilled since the chiral image of a L tetromino
is labeled as J, breaking label invariance. A 4-fold rotation, i.e.
$90^{\circ}$ rotation of the image can be expressed as a reordering of the
pixels described by the following index transformation:

$x^{\prime}_{i,j}=x_{n-1-j,i}$ (9)

A crucial step to build a label-invariant quantum model is to find the
representation of this operation acting on the Hilbert space of the quantum
circuit. To do so, first, the chosen data-encoding strategy will be described.

###  III.2 Angle encoding

Frequently named as qubit encoding, angle encoding [11, 48] is one of the most
simple data-encoding methods. It scales linearly in the number of required
qubits and can be implemented in constant time, making it one of the most fit
candidates to apply on NISQ devices. It consists on encoding the features
$x_{k}$ of the data point $x=(x_{0},\dotsi,x_{N})$ in the amplitudes of the
state-vector of each qubit:

$|x\rangle=\bigotimes_{k=0}^{N}\cos(x_{k})|0\rangle+\sin(x_{k})|1\rangle.$
(10)

This encoding can be performed by applying $R_{X}(x_{k})$ gates over each
qubit (see fig.3), where the pixel $x_{i,j}$ is mapped to $x_{k}$ vector value
with $k(i,j)=in+j$ index transformation.

Figure 3: Schematic example of the angle encoding method for a data-point with
4 features.

Due to the periodic nature of $R_{X}$ rotation gates, the classical features
ought to be scaled to fit in $x_{k}\in[-\pi,\pi]$. This method requires one
qubit per feature in the data-point, which can be reduced by a factor of two
encoding an additional feature in the complex phase of each qubit, following
the so-called dense angle encoding strategy. This improvement in the number of
required qubits $n_{q}$, which is still linear with respect to the number of
features $\mathcal{O}(N)$, is notably surpassed by the exponential reduction
achieved with amplitude encoding [49] $\mathcal{O}(\log_{2}(N))$.
Nevertheless, this improvement in the number of qubits comes along with an
exponential circuit depth and the impossibility of performing data re-
uploading. For this reason, the angle encoding method will be employed in this
work and a way to handle the issue of the number of qubits will be discussed
in subsequent stages.

###  III.3 Equivariant variational circuit

As a 4-fold rotation of an image can be performed by rearranging matrix
elements according to the index transformation described in eq.9 and every
matrix element is encoded in the quantum state of independent qubits, that
symmetry operation can be represented in terms of hermitian SWAP gates acting
on the Hilbert space. The mentioned index transformation can be achieved by
computing the transposed matrix of the original copy and reordering its
columns (see fig. 4).

(a) Transposition

(b) Column permutation

Figure 4: Rotation operation decomposed into a (4(a)) transposition and (4(b))
column permutation using SWAP operations. Example for an image of n=4.

More specifically, for the chosen data encoding method, the quantum unitary
that performs the action of the rotation symmetry element takes the following
form:

$\displaystyle U_{g}$ $\displaystyle=$
$\displaystyle\prod_{\begin{subarray}{c}i=1,\dotsi,n\\\
j=1,\dotsi,\lfloor\frac{n}{2}\rfloor\end{subarray}}SWAP(q_{k(i,j)},q_{k(i,n-j-1)})$
(11) $\displaystyle\prod_{\begin{subarray}{c}i=1,\dotsi,n\\\
j<i\end{subarray}}SWAP(q_{k(i,j)},q_{k(j,i)}).$

The unitary representation $U_{g}$ of the symmetry element $g\in G$ can be
used as a generator for the representations of the rest of the symmetry
elements in the group, operating by matrix multiplication with itself:

$\displaystyle U_{g^{2}}=U_{g}U_{g}=U_{g}^{2},$ $\displaystyle
U_{g^{3}}=U_{g}U_{g^{2}}=U_{g}^{3},$ $\displaystyle
U_{g^{4}}=U_{g}U_{g^{3}}=U_{g}^{4}.$ (12)

For this particular case, it is remarkable that $U_{g}^{3}$ represents a
rotation in the opposite direction and, coping with unitary operations,
$U_{g}^{3}=U_{g}^{-1}=U_{g}^{\dagger}$. Moreover, it is straightforward that
$U_{g}^{4}=\mathbb{I}_{2^{n_{q}}}$ represents the identity operation. Thus,
for this context, the symmetry group is the cyclic $C_{4}$ group and its
unitary representation is
$\mathcal{G}_{U}=\\{U_{g},U_{g}^{2},U_{g}^{\dagger},\mathbb{I}\\}$. Now it can
be verified that the equivariance condition described in eq. II.3 is satisfied
for the chosen encoding unitary.

In order to build the variational blocks $W_{l}(\theta_{l})$, parameterized
rotation gates will be used to rotate the quantum states of qubits in the
Bloch’s sphere and two-qubit CNOT gates to create entanglement between them.
For convenience, the generators of the quantum gates will be simmetrized,
which can be generated by Pauli operators, since
$[A,B]=0\Rightarrow[e^{A},B]=0$ :

$R_{X}(\theta)\equiv e^{-i\theta\frac{X}{2}}.$ (13)

Now, it is possible to symmetrize Pauli operators using the twirling formula.
If $X_{k}$ is defined as a Pauli gate $X$ applied on the qubit number $k$
(e.g. $X_{0}=\mathbb{I}\otimes\mathbb{I}\otimes\dotsi\otimes X$), the
following result is obtained:

$\displaystyle\mathcal{T}[X_{k}]$ $\displaystyle=$
$\displaystyle\frac{1}{4}(X_{0}+U_{g}X_{k}U_{g}^{\dagger}+U_{g}^{2}X_{k}{U_{g}^{\dagger}}^{2}+U_{g}^{3}X_{k}{U_{g}^{\dagger}}^{3})$
(14) $\displaystyle=$ $\displaystyle\frac{1}{4}(X_{k(i,j)}+X_{k(n-j-1,i)}+$
$\displaystyle X_{k(n-i-1,n-j-1)}+X_{k(j,n-i-1)}).$

Let orbit $\varphi$ denote the set of all different indices obtained under the
action of all symmetry operations for a given pair of indices $(i,j)$:

$\displaystyle\varphi(i,j)$ $\displaystyle=$
$\displaystyle\\{(i,j),(n-j-1,i),$ (15)
$\displaystyle(n-i-1,n-j-1),(j,n-i-1)\\}.$

Then, the symmetrization of a Pauli gate is equal for every same Pauli gate
applied on the qubits connected by the same orbit. This idea can be extended
to the parameterized rotation gates generated by those Pauli gates since
$[R_{X_{k}}(\theta),R_{X_{k^{\prime}}}(\theta^{\prime})]=0$:

$\displaystyle
e^{-i\frac{\theta}{8}(X_{\varphi_{0}(i,j)}+X_{\varphi_{1}(i,j)}+X_{\varphi_{2}(i,j)}+X_{\varphi_{3}(i,j)})}$
$\displaystyle=$ $\displaystyle
R_{X_{\varphi_{0}(i,j)}}(\theta^{\prime})R_{X_{\varphi_{1}(i,j)}}(\theta^{\prime})R_{X_{\varphi_{2}(i,j)}}(\theta^{\prime})R_{X_{\varphi_{3}(i,j)}}(\theta^{\prime}).$

Consequently, the same symmetrization is obtained for rotation gates with the
same rotation angle applied on qubits connected by the same orbit. This
symmetrization is equivalent for the Y and Z rotation axes. The same logic can
be followed for two-qubit gates, leading to this result: the same two-qubit
gate must be applied over the orbits of both qubits, which must be in
different orbits, to ensure the permutation condition of these gates with
respect to the representation of the symmetry group. To illustrate this, as a
reference to condensed matter physics, we will call unit cell to a set of
indices that don’t share orbits; this is, a set of points that are not related
by a symmetry operation (see fig. 5).

(a) Orbits

(b) Intra-cell

(c) Inter-cell

Figure 5:  A schematic image of different orbits (5(a)) and different unit
cells (5(b), 5(c)) represented by colors for an image of $n=4$. 5(b) and 5(c)
show the two-qubit gate architecture chosen to build an equivariant
variational unitary for an image of $n=4$. It is built out of intra-cell and
inter-cell $CNOT$ gates.

Apart from the constraints concluded from the symmetrization of two-qubit
gates, one has a certain freedom to choose their connectivity, i.e., on which
qubit it can apply them. Indeed, the choice of the unit cell shown in 5(b) and
5(c) is not unique. A circular connectivity scheme will be selected,which
entangles the qubits inside every unit cell. However, if different unit cells
are not connected between them, it will produce four independent quantum
circuits that would process just some phew pixels of the original image. To
avoid this scenario, some CNOTs can be added between unit cells and deal with
a more extensive set of correlated qubits.

Once a variational block has been built out of symmetrized gates, the sequence
of encoding and variational unitaries can be applied in an arbitrary number of
layers to build an equivariant data re-uploading circuit. It is of great
importance to make sure the chosen initial state of the quantum circuit is
invariant to the action of the symmetry group, as mentioned in sec. II.3. In
this case, it can verifed that $|0\rangle^{\otimes n_{q}}$ satisfies the
invariant condition. Finally, an invariant observable is the last ingredient
left to build a label-invariant model. A fit example for this problem can be
the sum of Pauli Z gates applied on qubits connected by an orbit:

$\displaystyle O$ $\displaystyle=$
$\displaystyle\frac{1}{4}(Z_{0}+Z_{3}+Z_{12}+Z_{15}),$ $\displaystyle
f(x,\theta)$ $\displaystyle=$
$\displaystyle\langle\psi(x,\theta)|O|\psi(x,\theta)\rangle.$ (17)

###  III.4 Simulation results

An equivariant quantum circuit has been trained to learn to classify tetromino
images described in III.1. This circuit employs parameterized quantum gates
symmetrized using the Twirling Formula, following the workflow detailed in
III.3. In consequence, a label-invariant model was obtained with respect to
the action of 4-fold rotations. To benchmark the power of the model and to
witness whether the geometric approach does provide any advantage in
comparison to conventional variational circuits, its performance was compared
to a widely used architecture called Basic Entangler Layers, based on [50]. To
make the comparison the fairest possible and since the latter architecture
makes use of more optimizable parameters per layer, which can induce an
improved approximation capacity, general rotation gates were used in the
Equivariant model, which can be decomposed into $3$ optimizable rotation
gates, up to a negligible global phase [51]:

$U=e^{i\alpha}R_{Z}(\beta)R_{Y}(\gamma)R_{Z}(\delta).$ (18)

It is true that these general rotation gates are more expressive than one-axis
rotation gates employed in the Basic Entangler model, and this could be a
possible source of a difference in performance. For this reason, a third model
type was compared, which we will call Non-equivariant model. Its quantum
circuit utilizes the same architecture as the Equivariant model’s but with
randomly generated orbits to break the symmetrization of gates and,
consequently, the equivariant condition. This way, we ensure we capture if any
improvement in the performance is due to the choice of gates or if it is,
indeed, a consequence of the geometric approach. See tab. 1 to compare the
number of optimizable parameters in each model.

Model | Parameters
---|---
Equivariant | $3n_{\varphi}n_{l}=12n_{l}$
Non-equivariant | $3n_{\varphi}n_{l}=12n_{l}$
Basic Entangler | $n_{q}n_{l}=16n_{l}$
Table 1: The number of optimizable parameters in the quantum circuits.
$n_{\varphi}$ corresponds to the number of different orbits in an image,
$n_{l}$ is the number of data re-uploading layers, and $n_{q}$ corresponds to
the number of qubits in the circuit. For an image of $n=4$, $n_{\varphi}=4$
and $n_{q}=16$. Figure 6: Loss landscape of three model types along the
training process where $10$ models were trained for each type. The shaded
areas are spread between the maximum and minimum costs of each model type, and
the dashed lines represent the mean cost of each model type. All models apply
$n_{l}=5$ re-uploading layers.

All three models were trained in a noiseless environment using the default
simulator provided by Pennylane [52], and the parameter update is performed
using Adam Optimizer [53]. In fig. 6, it can observed that the Equivariant
model reaches the lowest values in the loss landscape for the same number of
iterations. Even if the Basic Entangler performs a steeper descent in early
steps, it reaches a plateau where the cost nearly stops decreasing. The other
two models, however, keep reducing the cost until the last optimization step,
although their descent tends to flatten in the last iterations.

The loss function, defined as the mean quadratic distance between the
circuit’s predictions and the true labels,
$\mathcal{L}(x,\hat{y},\theta)=\frac{1}{N}\sum_{i}^{N}(f(x^{(i)},\theta)-\hat{y}^{(i)})^{2}$,
does not capture accurately the performance of the models, and lower costs do
not always result in better predictions. To evaluate and compare the
predictive power of each model type, different metrics are used, such as the
accuracy, precision, recall and f1 score [54]. In fig. 7, the performance of
various trained models is compared. It is remarkable how the equivariant
models outperform the other two model types both in fitting the data and in
generalization.

Figure 7: Violin plots of F1 score metrics of trained models. Ten different
models of each type were trained with the same hyperparameters ($n_{l}=5$),
and the performance was evaluated over the training and test set.

Although the results shown in fig.7 are favorable for the geometric solution,
the performances are overall improvable, especially in their generalization
capability. A classical procedure to handle this matter could be augmenting
the data set size; nonetheless, the function classes the models can generate
will be restricted by the number of data re-uploading layers applied in the
quantum circuits. In consequence, the models will reach a performance ceiling
that couldn’t be improved even with more training steps. Performing additional
data re-uploading layers will give the quantum circuit the capacity to
generate higher-order trigonometrical polynomials of the input data features
and, in conclusion, improve its approximation power. In Tab.2, performance
metrics for different numbers of data re-uploading layers are shown.

Model Type | $n_{l}$ | F1 | Accuracy | Precision | Recall
---|---|---|---|---|---
Equiv | 2 | 0.59/0.47 | 0.66/0.52 | 0.72/0.57 | 0.51/0.40
Non-equiv | 2 | 0.53/0.50 | 0.61/0.53 | 0.65/0.58 | 0.45/0.44
Basic Entangler | 2 | 0.51/0.47 | 0.55/0.47 | 0.54/0.50 | 0.49/0.45
Equiv | 4 | 0.78/0.63 | 0.81/0.64 | 0.85/0.68 | 0.73/0.59
Non-equiv | 4 | 0.55/0.41 | 0.62/0.48 | 0.66/0.52 | 0.47/0.34
Basic Entangler | 4 | 0.47/0.37 | 0.55/0.46 | 0.55/0.48 | 0.41/0.31
Equiv | 6 | 0.79/0.67 | 0.82/0.70 | 0.87/0.77 | 0.73/0.61
Non-equiv | 6 | 0.63/0.45 | 0.68/0.54 | 0.71/0.60 | 0.57/0.37
Basic Entangler | 6 | 0.67/0.46 | 0.73/0.52 | 0.80/0.56 | 0.59/0.41
Equiv | 8 | 0.69/0.48 | 0.75/0.58 | 0.85/0.72 | 0.58/0.36
Non-equiv | 8 | 0.63/0.50 | 0.68/0.57 | 0.75/0.69 | 0.54/0.41
Basic Entangler | 8 | 0.66/0.47 | 0.72/0.52 | 0.79/0.58 | 0.56/0.41
Equiv | 10 | 0.98/0.93 | 0.98/0.93 | 1.00/1.00 | 0.96/0.88
Non-equiv | 10 | 0.73/0.54 | 0.78/0.60 | 0.82/0.65 | 0.67/0.47
Basic Entangler | 10 | 0.87/0.73 | 0.87/0.74 | 0.89/0.81 | 0.85/0.67

Table 2: Performance metrics for each model type and per different number of
data re-uploading layers. Each value in the table represents the mean
performance metric of 5 different models, evaluated over the training set
(left) and the test set (right).

The results in Tab.2 show, on one hand, that, in general, the Equivariant
model achieves better performances both in the training and in the test sets.
Moreover, it is noticeable how increasing the number of data re-uploading
layers tends to improve the approximation ability of all three models; this is
well illustrated in Fig.8.

Figure 8: Mean $F1$ scores over the test set of the three model types (10
modes per type) in function of the number of data re-uploading layers $n_{l}.$

###  III.5 Scaling the problem size

The data set employed in the previous sections represents a toy case created
for the comparison of the capabilities of the geometric approach with respect
to the traditional variational quantum models. The resolution of images in
real computer vision problems, however, is far from the phew pixels the
current NISQ-era devices can handle due to their limited circuit size. A way
to address this issue is to use a another encoding strategy, different from
the one proposed in this work, that optimizes the number of features encoded
per qubit. In [44, 43], for example, Amplitude Encoding is used the encode the
images into the quantum circuits. Another way to cope with this matter, rather
than waiting to larger-scale quantum devices to be built and implement the
proposed architecture for larger images, is to reduce the dimensionality of
the input.

There are several feature reduction techniques such as PCA, as the main lineal
transformation method, or autoencoders among the non-linear transformations
[55]. In this work, we propose to use trainable convolutional layers to
extract features in the latent space, preserving the spatial arrangement of
the features in the original image. Nevertheless, if the model is wanted to be
label-invariant, the classical convolution layers must be equivariant to the
operations in the symmetry group. A two-dimensional discrete convolution
between an image $x$ and a filter $w$ can be described by the following
equation, assuming square images and filters [56]:

$(x\ast w)[i,j]=\sum_{s=-r}^{r}\sum_{t=-r}^{r}x[s,t]w[i-s,j-t],$ (19)

where $r=(n_{f}-1)/2$ for odd $n_{f}$ and $r=n_{f}/2$ for even $n_{f}$, being
$n_{f}$ the size of the filter. This expression can be extended to consider
RGB channels, but in general, the operation is not equivariant to image
rotations described by eq. 9. For this reason, following the approach of [57],
a rotation-equivariant convolution operation is defined by computing a
convolution between an image and all rotated versions of the original filter,
and pooling all results (see fig. 9).

Figure 9: Rotation equivariant convolution operation.

These equivariant convolutional layers can be stacked arbitrarily in $n_{c}$
channels, each one of them with a different filter so that they learn to focus
on specific characteristics of the images. This way, if $n_{x}$ is the size
(number of pixels) of the rows of the square image, $n_{w}$ the size of the
filter, $s$ the stride and $p$ the padding, the size of the convoluted image
will be [58]:

$n^{\prime}_{x}=\frac{n_{x}+2p-n_{w}}{s}+1$ (20)

For this case, $p=0$ and $s=1$ will be taken, so
$n^{\prime}_{x}=n_{x}-n_{w}+1$. Taking the data set as a tensor of $N$ data-
points allows us to apply the equivariant convolution operation over the whole
batch, performing the following dimensional reduction:

Input tensor: $\displaystyle(N,n_{x},n_{x},n_{c})$ Filter:
$\displaystyle(n_{w},n_{w},n_{c},n^{\prime}_{c})$ Output tensor
$\displaystyle(N,n^{\prime}_{x},n^{\prime}_{x},n^{\prime}_{c}).$ (21)

Accordingly, the number of channels of each convolutional layer and the size
of the filters can be tuned to reduce the number of features and adapt the
input to the size of the available quantum circuit. In this case, to connect
the classical part to the equivariant quantum classifier described in the
previous sections, the original images will be reduced to $4\times 4$, single
channel square images. It is postulated that while the number of qubits keeps
increasing and the connectivity between them condensing, the size of the
images that can be embedded directly into the quantum classifier will increase
while decreasing the dependency on the classical convolutional part.

This algorithm was tested on two public datasets. The first one, Fashion
MNIST111https://huggingface.co/datasets/fashion_mnist, contains $28\times 28$
grayscale images. For binary classification, labels $0$ (T-shirts) and $1$
(trousers) have been taken and for a dataset of 200 images, convolution
filters of dimensions of $n_{w}=\\{11,11,3,3\\}$ and $n_{l}=2$ data re-
uploading layers, a mean performance of $F1_{\text{train}}=0.82\pm 0.09$ and
$F1_{\text{test}}=0.79\pm 0.11$ was obtained. The second dataset, 11k
Hands222https://sites.google.com/view/11khands [59], has been chosen as an
illustrative example of a rotational symmetric problem where reflection
symmetry is not present. It contains $1600\times 1200$ RGB images of hands
that have been reduced to a resolution of $64\times 64$ before applying the
hybrid algorithm. 500 images of left and right hands have been taken for
binary classification, and for filters of dimensions
$n_{w}=\\{31,7,11,11,3,3\\}$, a dataset of 500 images and $n_{l}=2$ re-
uploading layers, a perfect classification was achieved in all the cases, both
in the training and testing sets.

##  IV Conclusion

This work has shown a methodology to build a label-invariant VQC for image
classification with $C_{4}$ rotation label-symmetry. This has been achieved by
symmetrizing the quantum gates using the Twirling formula and choosing the
right initial state and observable. The model has been benchmarked in a
synthetic data set, and compared against two quantum models with different
architectures. The first one is a standard architecture, widely used in the
literature; the second uses the same gates, initial state, and observable as
the label-invariant model but with randomized orbits, breaking any symmetry-
preserving relation.

Experimental results in simulators have shown that the equivariant quantum
circuit obtains better performances than the other two models, coinciding with
the expected results from the literature of GQML. It reaches lower minima than
the two other models and the performance metrics, both in the training and
test sets, are higher on average. The difference in the performance between
the Equivariant and Non-equivariant models reveals that the improvement in
classification comes from the geometric approach since the same gates,
observable, and number of parameters have been used, and the only change is
the qubits they are acting on. Alternatively, the importance of data re-
uploading has been observed in the approximation capabilities of every model,
increasing the performance metrics with the number of re-uploading layers.

The quantum models were designed for a synthetic toy data set of phew
grayscale pixels, considering the reduced number of available qubits and low
connectivity in NISQ devices, and the computational cost of simulating them.
An equivariant convolution operation has been proposed to implement the
quantum circuits for higher-resolution images with color channels while
preserving the label-invariance of the model. This way, the weights of
convolution filters have been trained to learn the most essential features of
the image, which are introduced as input to the quantum circuit. This hybrid
algorithm has been tested on two public datasets to test its feasibility.

Taking into account that all the experiments have been made in an ideal,
noise-free simulator, we consider it of great interest to execute the models
in real quantum hardware and observe their behavior with noise. In addition,
the chosen data encoding method embeds a single grayscale value per qubit. The
most straightforward solution to embed color channels could be to execute
parallel equivariant circuits, one for each channel, entangled between them.
Another approach could be to use the imaginary phase of qubits and map every
point of Bloch’s sphere to a two-dimensional colormap using dense angle
encoding. We delegate to future investigations to research qubit-efficient
data encoding methods that consider color channels and are equivariant to
image rotations. The choice of the encoding will determine the unitary
representation of the symmetry group and, hence, the symmetrization of gates.
Extending the model for multi-class classification could be another task for
future work, choosing the right observable that returns multiple labels and
preserves the symmetry relations of the data.

Finally, we believe variational quantum algorithms have great potential as an
early application of NISQ devices, although the anticipated quantum advantage
may not come in the shape and flavor it is expected [60]. All in all, lots of
work is to be done, and the geometric framework seems a powerful tool to
design new architectures for machine learning models, even beyond quantum
algorithms. We hope this work serves as an example of the methodology to build
variational quantum classifiers with geometric inductive bias and encourages
future research to propose new and innovative applications of quantum
technologies.

Acknowledgements.- We acknowledge Ikerlan, DIPC, Ikerbasque, Basque
Government, Diputación de Gipuzkoa and EIC for support, as well as discussions
with the teams from Multiverse Computing and Ikerlan.

## References

  * Fedorov  _et al._ [2022] A. K. Fedorov, N. Gisin, S. M. Beloussov, and A. I. Lvovsky, “Quantum computing at the quantum advantage threshold: a down-to-business review,”  (2022), arXiv:2203.17181 [physics, physics:quant-ph].
  * Abbas  _et al._ [2023] A. Abbas, A. Ambainis, B. Augustino, A. BÃ¤rtschi, H. Buhrman, C. Coffrin, G. Cortiana, V. Dunjko, D. J. Egger, B. G. Elmegreen, N. Franco, F. Fratini, B. Fuller, J. Gacon, C. Gonciulea, S. Gribling, S. Gupta, S. Hadfield, R. Heese, G. Kircher, T. Kleinert, T. Koch, G. Korpas, S. Lenk, J. Marecek, V. Markov, G. Mazzola, S. Mensa, N. Mohseni, G. Nannicini, C. O’Meara, E. P. Tapia, S. Pokutta, M. Proissl, P. Rebentrost, E. Sahin, B. C. B. Symons, S. Tornow, V. Valls, S. Woerner, M. L. Wolf-Bauwens, J. Yard, S. Yarkoni, D. Zechiel, S. Zhuk, and C. Zoufal, “Quantum Optimization: Potential, Challenges, and the Path Forward,”  (2023), arXiv:2312.02279 [quant-ph].
  * Au-Yeung  _et al._ [2023] R. Au-Yeung, B. Camino, O. Rathore, and V. Kendon, “Quantum algorithms for scientific applications,”  (2023), arXiv:2312.14904 [physics, physics:quant-ph].
  * Orús  _et al._ [2019] R. Orús, S. Mugel, and E. Lizaso, Reviews in Physics 4, 100028 (2019).
  * Aizpurua  _et al._ [2023] B. Aizpurua, P. Bermejo, J. E. Martinez, and R. Orús, “Hacking cryptographic protocols with advanced variational quantum attacks,” (2023), arXiv:2311.02986 [quant-ph] .
  * Bermejo and Orús [2023] P. Bermejo and R. Orús, Scientific Reports 13, 9840 (2023).
  * Bermejo and Orús [2022] P. Bermejo and R. Orús, “Variational quantum continuous optimization: a cornerstone of quantum mathematical analysis,” (2022), arXiv:2210.03136 [quant-ph] .
  * Bermejo  _et al._ [2023] P. Bermejo, B. Aizpurua, and R. Orús, “Improving gradient methods via coordinate transformations: Applications to quantum machine learning,” (2023), arXiv:2304.06768 [quant-ph] .
  * Wittek [2014] P. Wittek,  _Quantum machine learning: what quantum computing means todata mining_ , 1st ed. (Elsevier, AP, Amsterdam, 2014).
  * Biamonte  _et al._ [2017] J. Biamonte, P. Wittek, N. Pancotti, P. Rebentrost, N. Wiebe, and S. Lloyd, Nature 549, 195 (2017), arXiv:1611.09347 [cond-mat, physics:quant-ph, stat].
  * Schuld and Petruccione [2021] M. Schuld and F. Petruccione,  _Machine Learning with Quantum Computers_, Quantum Science and Technology (Springer International Publishing, Cham, 2021).
  * Liu  _et al._ [2023a] J. Liu, F. Wilde, A. A. Mele, L. Jiang, and J. Eisert, “Stochastic noise can be helpful for variational quantum algorithms,”  (2023a), arXiv:2210.06723 [quant-ph].
  * Wang and Liu [2024] Y. Wang and J. Liu, “Quantum Machine Learning: from NISQ to Fault Tolerance,”  (2024), arXiv:2401.11351 [quant-ph, stat].
  * Cerezo  _et al._ [2021a] M. Cerezo, A. Arrasmith, R. Babbush, S. C. Benjamin, S. Endo, K. Fujii, J. R. McClean, K. Mitarai, X. Yuan, L. Cincio, and P. J. Coles, Nature Reviews Physics 3, 625 (2021a), arXiv:2012.09265 [quant-ph, stat].
  * Havlicek  _et al._ [2019] V. Havlicek, A. D. Córcoles, K. Temme, A. W. Harrow, A. Kandala, J. M. Chow, and J. M. Gambetta, Nature 567, 209 (2019), arXiv:1804.11326 [quant-ph, stat].
  * Farhi and Neven [2018] E. Farhi and H. Neven, “Classification with Quantum Neural Networks on Near Term Processors,”  (2018), arXiv:1802.06002 [quant-ph].
  * Rodriguez-Grasa  _et al._ [2024] P. Rodriguez-Grasa, Y. Ban, and M. Sanz, “Training embedding quantum kernels with data re-uploading quantum neural networks,”  (2024), arXiv:2401.04642 [quant-ph].
  * Pérez-Salinas  _et al._ [2019] A. Pérez-Salinas, A. Cervera-Lierta, E. Gil-Fuster, and J. I. Latorre, Quantum 4 (2019), 10.22331/q-2020-02-06-226, arXiv: 1907.02085v3 Publisher: Verein zur Forderung des Open Access Publizierens in den Quantenwissenschaften.
  * Schuld  _et al._ [2020a] M. Schuld, R. Sweke, and J. J. Meyer, Physical Review A 103 (2020a), 10.1103/PhysRevA.103.032430, arXiv: 2008.08605v2 Publisher: American Physical Society.
  * McClean  _et al._ [2018] J. R. McClean, S. Boixo, V. N. Smelyanskiy, R. Babbush, and H. Neven, Nature Communications 9, 4812 (2018), number: 1 Publisher: Nature Publishing Group.
  * Holmes  _et al._ [2022] Z. Holmes, K. Sharma, M. Cerezo, and P. J. Coles, PRX Quantum 3, 010313 (2022).
  * Larocca  _et al._ [2022a] M. Larocca, P. Czarnik, K. Sharma, G. Muraleedharan, P. J. Coles, and M. Cerezo, Quantum 6, 824 (2022a), publisher: Verein zur Förderung des Open Access Publizierens in den Quantenwissenschaften.
  * Larocca  _et al._ [2023] M. Larocca, N. Ju, D. García-Martín, P. J. Coles, and M. Cerezo, Nature Computational Science 3, 542 (2023), number: 6 Publisher: Nature Publishing Group.
  * Cerezo  _et al._ [2021b] M. Cerezo, A. Sone, T. Volkoff, L. Cincio, and P. J. Coles, Nature Communications 12, 1791 (2021b), number: 1 Publisher: Nature Publishing Group.
  * Ragone  _et al._ [2023a] M. Ragone, B. N. Bakalov, F. Sauvage, A. F. Kemper, C. O. Marrero, M. Larocca, and M. Cerezo, “A Unified Theory of Barren Plateaus for Deep Parametrized Quantum Circuits,”  (2023a), arXiv:2309.09342 [quant-ph].
  * Fontana  _et al._ [2023] E. Fontana, D. Herman, S. Chakrabarti, N. Kumar, R. Yalovetzky, J. Heredge, S. H. Sureshbabu, and M. Pistoia, “The Adjoint Is All You Need: Characterizing Barren Plateaus in Quantum Ans\”atze,”  (2023), arXiv:2309.07902 [quant-ph].
  * Park and Killoran [2024] C.-Y. Park and N. Killoran, Quantum 8, 1239 (2024).
  * Liu  _et al._ [2023b] H.-Y. Liu, T.-P. Sun, Y.-C. Wu, Y.-J. Han, and G.-P. Guo, New Journal of Physics 25, 013039 (2023b).
  * Sack  _et al._ [2022] S. H. Sack, R. A. Medina, A. A. Michailidis, R. Kueng, and M. Serbyn, PRX Quantum 3, 020365 (2022), publisher: American Physical Society.
  * Nguyen  _et al._ [2022] Q. T. Nguyen, L. Schatzki, P. Braccia, M. Ragone, P. J. Coles, F. Sauvage, M. Larocca, and M. Cerezo, “Theory for Equivariant Quantum Neural Networks,”  (2022).
  * Ragone  _et al._ [2023b] M. Ragone, P. Braccia, Q. T. Nguyen, L. Schatzki, P. J. Coles, F. Sauvage, M. Larocca, and M. Cerezo, “Representation Theory for Geometric Quantum Machine Learning,”  (2023b), arXiv:2210.07980 [quant-ph, stat].
  * Meyer  _et al._ [2023] J. J. Meyer, M. Mularski, E. Gil-Fuster, A. A. Mele, F. Arzani, A. Wilms, and J. Eisert, PRX Quantum 4, 010328 (2023), arXiv:2205.06217 [quant-ph].
  * Larocca  _et al._ [2022b] M. Larocca, F. Sauvage, F. M. Sbahi, G. Verdon, P. J. Coles, and M. Cerezo, PRX Quantum 3, 030341 (2022b).
  * Schatzki  _et al._ [2022] L. Schatzki, M. Larocca, Q. T. Nguyen, F. Sauvage, and M. Cerezo, “Theoretical Guarantees for Permutation-Equivariant Quantum Neural Networks,”  (2022).
  * Bronstein  _et al._ [2021] M. M. Bronstein, Joan Bruna, T. Cohen, and P. Veličković, “Geometric Deep Learning: Grids, Groups, Graphs, Geodesics, and Gauges,”  (2021), arXiv:2104.13478 [cs, stat].
  * Bowles  _et al._ [2023] J. Bowles, V. J. Wright, M. Farkas, N. Killoran, and M. Schuld, “Contextuality and inductive bias in quantum machine learning,”  (2023), arXiv:2302.01365 [quant-ph].
  * Tüysüz  _et al._ [2024] C. Tüysüz, S. Y. Chang, M. Demidik, K. Jansen, S. Vallecorsa, and M. Grossi, “Symmetry breaking in geometric quantum machine learning in the presence of noise,”  (2024), arXiv:2401.10293 [quant-ph].
  * Zheng  _et al._ [2022] H. Zheng, G. S. Ravi, H. Wang, K. Setia, F. T. Chong, and J. Liu,  (2022), arXiv: 2211.12711.
  * West  _et al._ [2023a] M. T. West, J. Heredge, M. Sevior, and M. Usman, “Provably Trainable Rotationally Equivariant Quantum Machine Learning,”  (2023a), arXiv:2311.05873 [quant-ph].
  * East  _et al._ [2023] R. D. P. East, G. Alonso-Linaje, and C.-Y. Park, “All you need is spin: SU(2) equivariant variational quantum circuits based on spin networks,”  (2023), arXiv:2309.07250 [cond-mat, physics:quant-ph, stat].
  * Heredge  _et al._ [2023] J. Heredge, C. Hill, L. Hollenberg, and M. Sevior,  (2023), arXiv: 2304.03601.
  * Skolik  _et al._ [2022] A. Skolik, M. Cattelan, S. Yarkoni, T. Bäck, and V. Dunjko,  (2022), 10.1038/s41534-023-00710-y, arXiv: 2205.06109.
  * West  _et al._ [2023b] M. T. West, M. Sevior, and M. Usman, Machine Learning: Science and Technology 4, 035027 (2023b), arXiv:2212.00264 [quant-ph].
  * Chang  _et al._ [2023] S. Y. Chang, M. Grossi, B. L. Saux, and S. Vallecorsa, “Approximately Equivariant Quantum Neural Network for $p4m$ Group Symmetries in Images,”  (2023), arXiv:2310.02323 [quant-ph].
  * Vapnik [1991] V. Vapnik, in  _Advances in Neural Information Processing Systems_, Vol. 4 (Morgan-Kaufmann, 1991).
  * Wierichs  _et al._ [2022] D. Wierichs, J. Izaac, C. Wang, and C. Y.-Y. Lin, Quantum 6, 677 (2022), publisher: Verein zur Förderung des Open Access Publizierens in den Quantenwissenschaften.
  * Davvaz [2021] B. Davvaz,  _Groups and Symmetry: Theory and Applications_ (2021).
  * [48] R. Larose and B. Coyle, ArXiv: 2003.01695v1.
  * Mottonen  _et al._ [2004] M. Mottonen, J. J. Vartiainen, V. Bergholm, and M. M. Salomaa, “Transformation of quantum states using uniformly controlled rotations,”  (2004), arXiv:quant-ph/0407010.
  * Schuld  _et al._ [2020b] M. Schuld, A. Bocharov, K. Svore, and N. Wiebe, Physical Review A 101, 032308 (2020b), arXiv:1804.00633 [quant-ph].
  * Nielsen and Chuang [2010] M. A. Nielsen and I. L. Chuang, Quantum Computation and Quantum Information (2010), 10.1017/CBO9780511976667, publisher: Cambridge University Press ISBN: 9780511976667.
  * Bergholm  _et al._ [2018] V. Bergholm, J. Izaac, M. Schuld, C. Gogolin, S. Ahmed, V. Ajith, M. S. Alam, G. Alonso-Linaje, B. AkashNarayanan, A. Asadi, J. M. Arrazola, U. Azad, S. Banning, C. Blank, T. R. Bromley, B. A. Cordier, J. Ceroni, A. Delgado, O. Di Matteo, A. Dusko, T. Garg, D. Guala, A. Hayes, R. Hill, A. Ijaz, T. Isacsson, D. Ittah, S. Jahangiri, P. Jain, E. Jiang, A. Khandelwal, K. Kottmann, R. A. Lang, C. Lee, T. Loke, A. Lowe, K. McKiernan, J. J. Meyer, J. A. Montañez-Barrera, R. Moyard, Z. Niu, L. J. O’Riordan, S. Oud, A. Panigrahi, C.-Y. Park, D. Polatajko, N. Quesada, C. Roberts, N. Sá, I. Schoch, B. Shi, S. Shu, S. Sim, A. Singh, I. Strandberg, J. Soni, A. Száva, S. Thabet, R. A. Vargas-Hernández, T. Vincent, N. Vitucci, M. Weber, D. Wierichs, R. Wiersema, M. Willmann, V. Wong, S. Zhang, and N. Killoran,  (2018), arXiv: 1811.04968.
  * Kingma and Ba [2017] D. P. Kingma and J. Ba, “Adam: A Method for Stochastic Optimization,”  (2017), arXiv:1412.6980 [cs].
  * Tharwat [2020] A. Tharwat, Applied Computing and Informatics 17, 168 (2020), publisher: Emerald Publishing Limited.
  * van der Maaten  _et al._ [2007] L. van der Maaten, E. Postma, and H. Herik, Journal of Machine Learning Research - JMLR 10 (2007).
  * Gonzalez and Woods [2018] R. C. Gonzalez and R. E. Woods,  _Digital image processing_ (Pearson, New York, NY, 2018).
  * Cohen and Welling [2016] T. S. Cohen and M. Welling, “Group Equivariant Convolutional Networks,”  (2016), arXiv:1602.07576 [cs, stat].
  * Dumoulin and Visin [2018] V. Dumoulin and F. Visin, “A guide to convolution arithmetic for deep learning,”  (2018), arXiv:1603.07285 [cs, stat].
  * Afifi [2019] M. Afifi, Multimedia Tools and Applications (2019), 10.1007/s11042-019-7424-8.
  * Cerezo  _et al._ [2023] M. Cerezo, M. Larocca, D. García-Martín, N. L. Diaz, P. Braccia, E. Fontana, M. S. Rudolph, P. Bermejo, A. Ijaz, S. Thanasilp, E. R. Anschuetz, and Z. Holmes, “Does provable absence of barren plateaus imply classical simulability? Or, why we need to rethink variational quantum computing,”  (2023), arXiv:2312.09121 [quant-ph, stat].


# Optical Identification of the Shortest-Period Spider Pulsar System M71E

Zhuokai Liu  Department of Astronomy, School of Physics, Peking University,
Yiheyuan Rd. 5, Haidian District, Beijing, China, 100871
Kavli Institute of Astronomy and Astrophysics, Peking University, Yiheyuan Rd.
5, Haidian District, Beijing, China, 100871
Subo Dong  Department of Astronomy, School of Physics, Peking University,
Yiheyuan Rd. 5, Haidian District, Beijing, China, 100871
Kavli Institute of Astronomy and Astrophysics, Peking University, Yiheyuan Rd.
5, Haidian District, Beijing, China, 100871

###### Abstract

M71E is a spider pulsar (i.e., a millisecond pulsar with a tight binary
companion) with the shortest known orbital period of $P=53.3$ min discovered
by Pan et al. (2023). Their favored evolutionary model suggests that it
bridges between two types of spider pulsars, namely, it descended from a
“redback” and will become a “black widow”. Using Hubble Space Telescope (HST)
archival imaging data, we report the first optical identification of its
companion COM-M71E. The HST and pulsar timing coordinates are in excellent
agreement (within $\sim 10$ mas). If M71E is associated with the globular
cluster M71, our measured brightness of COM-M71E ($m_{\rm F606W}\approx 25.3$)
is broadly consistent with the expectation from Pan et al. (2023)’s preferred
binary evolutionary model of a stripped dwarf companion, while it is also
compatible with an ultra-low-mass degenerate companion. Future multi-
wavelength photometric and spectroscopic observations can characterize the
companion and test the evolutionary scenarios.

Millisecond pulsars (1062) — Binary pulsars (153) — Optical identification
(1167) — Globular clusters (656)

††software: hst1pass, MAAT

##  1 Introduction

Millisecond pulsars (MSPs), rapidly-spinning pulsars (PSRs) found in either
the Galactic field or globular clusters (GCs), are thought to be “recycled”
old neutron stars, which have spun up via accreting mass from binary
companions (Alpar et al., 1982; Bhattacharya & van den Heuvel, 1991). In the
classic picture, the accretion process occurs as low-mass X-ray binaries
(LMXBs), after which MSPs are left with companions evolving into He white
dwarfs (WDs). Most binary MSPs have WD companions following the mass-period
relation expected from binary evolution models leading to He WDs (Tauris &
Savonije, 1999). However, the spider pulsars, with short orbital periods
($P\lesssim 1$ d) and primarily found in eclipsing systems, are exceptions to
such a relation (see, e.g., Roberts, 2013): The “black widows” (BWs) have very
low-mass companions ($M_{\rm c}\ll 0.1\,M_{\odot}$), and in contrast, the
companions of “redbacks” (RBs) are comparatively more massive ($M_{\rm c}\sim
0.1-0.4\,M_{\odot}$). The companion’s evolution is affected by
irradiation/ablation by the pulsar wind – BWs may eventually become isolated
MSPs after the evaporating companions are obliterated, and Benvenuto et al.
(2014) suggest that BWs descend from ablated RBs. The formation mechanisms of
spider pulsars are debated (e.g., King et al., 2003, 2005; Chen et al., 2013;
Benvenuto et al., 2014; Jia & Li, 2015; Ginzburg & Quataert, 2020).

By analyzing the pulsar timing data from the Five-hundred-meter Aperture
Spherical radio Telescope (FAST), Pan et al. (2023) (hereafter Pan23) find
that M71E (a.k.a., PSR J1953+1844, discovered by Han et al. 2021 in the field
of globular cluster M71) is an extraordinary non-eclipsing spider pulsar,
which has a record-breakingly short orbital period of $P=53.3$ min and a very
small mass function of $2.3\times 10^{-7}M_{\odot}$. Pan23’s analysis favors a
stripped dwarf companion with $M_{c}=0.047-0.097\,M_{\odot}$, estimated by
imposing mass-radius relations for brown dwarfs (BDs)/low-mass dwarfs (Burrows
et al., 1993) and restricting its radius within the Roche lobe radius
(Eggleton, 1983). Such a moderate mass and its very short period suggest that
it evolved from a RB and will evolve into a BW. Their favored binary evolution
model based on Chen et al. (2013) expects that the companion has an effective
temperature $T_{\rm eff}\sim 4500$ K and luminosity $L\sim 3.25\times
10^{-3}\,L_{\odot}$, which is much fainter than their upper limit on the
optical brightness ($>22$ mag) from non detections on SDSS images.

In this paper, we detect the optical companion (hereafter COM-M71E) of M71E on
an archival HST image. Our detection adds to a handful of optically identified
GC spiders with very low-mass companions (e.g., Pallanca et al., 2014;
Cadelano et al., 2015).

Figure 1: A $5\arcsec\times 5\arcsec$ section of the HST ACS F606W image
around COM-M71E. The red circle is centered at the pulsar timing position of
M71E, which agrees with the HST F606W position within $\sim 10$ mas, and the
circle’s radius is 10 times of the HST astrometric error.

##  2 Optical Observations & Data Analysis

We query the Mikulski Archive for Space Telescopes for archival HST images
covering the location of M71E, which is $\sim 2.5\arcmin$ away from the center
of M71. We find a single image (see Figure 1 for a section centered on the
location of M71E derived from pulsar timing) taken with the F606W filter and
an exposure time of 339 s by the Advanced Camera for Surveys (ACS) on UT
November 22, 2021 (mid exposure time at MJD = 59540.11878) associated with
program 16871111All the HST data used in this paper can be found in MAST:
http://dx.doi.org/10.17909/xb70-ee36 (catalog 10.17909/xb70-ee36) (PI: J.
Anderson; Anderson et al. 2021), which was proposed to make astrometric
correction of the charge transfer efficiency (CTE) effects. Visual inspections
near the M71E position reveal a faint point source that is $\sim 2\arcsec$
away from a bright star.

We employ the software hst1pass222https://www.stsci.edu/~jayander/HST1PASS/
(Anderson, 2022) based on the effective point-spread function (ePSF) method
(Anderson & King, 2000) to perform astrometric and photometric analysis using
the pixel-by-pixel CTE-corrected data product (i.e., the image with “_flc”
suffix). We set the parameters ${\rm FMIN}$ = 100 and $\rm HMIN=3$ to allow
detecting faint sources and finding stars near bright neighbors, respectively.
The faint source close to M71E visible to by-eye inspections is detected by
hst1pass. To investigate possible issues of extracting a faint source in the
proximity of a bright star, we inject 10 artificial stars with the same
instrumental magnitude as the source of interest in the vicinity of a bright
star sharing a similar background level and gradient, and the injections are
at random sub-pixel positions. We recover the astrometric position within
$\sim 5$ mas and magnitude within $\sim 0.1$ mag, suggesting that the source
extraction is reliable.

We subsequently compare the resulting star catalog with the HST ACS F606W
results of M71 taken on UT May 12, 2006 in Sarajedini et al. (2007). We
perform zero-point magnitude calibration using the common stars, and we
determine that COM-M71E has $m_{\rm F606W}=25.34\pm 0.31$ in the VEGAmag
system. We also estimate the relative astrometric uncertainties as a function
of magnitude from epoch-to-epoch comparisons. Then we perform astrometric
calibration using common stars with the Gaia Data Release 3 (DR3) catalog
(Gaia Collaboration et al., 2016, 2023). We apply 2nd-order 2D-polynomial
transformations of the distortion-corrected “(r,d)” coordinates from hst1pass
to the Gaia DR3 frame and translate the Gaia coordinates to the HST epoch by
considering the proper motions. Our derived J2000.0 position of COM-M71E is
(RA, Dec) $=(19^{\rm h}53^{\rm m}37.947^{\rm s}\pm 0.001^{\rm s},18\arcdeg
44\arcmin 54.32\arcsec\pm 0.01\arcsec)$, and the uncertainties are estimated
by combining the relative error derived from HST epoch-to-epoch comparisons
and the Gaia calibration error. The comparison of the HST position with
Pan23’s pulsar timing position at (RA, Dec) $=(19^{\rm h}53^{\rm
m}37.9464^{\rm s}\pm 0.0001^{\rm s},18\arcdeg 44\arcmin 54.310\arcsec\pm
0.002\arcsec)$ shows an excellent agreement (within $\sim 10$ mas) at the
$<1\,\sigma$ level. Therefore, the HST detection of COM-M71E is a secure
optical identification.

We assume that M71E is associated with M71, which has a distance of $D=4.0$
kpc (Baumgardt & Vasiliev, 2021). We adopt an extinction correction of $A_{\rm
F606W}=0.6$ estimated using the M71 reddening value from Dotter et al. (2010)
and the extinction coefficient from Sirianni et al. (2005). Then we obtain
COM-M71’s absolute magnitude $M_{\rm F606W}=11.7\pm 0.3$.

Figure 2: Comparison between Pan23’s favored binary evolutionary model and the
observational constraints from HST photometry of COM-M71E in the Hertzsprung-
Russell (HR) diagram. The blue curve is the evolutionary track presented in
Panel (A) of Extended Data Figure 4 of Pan23. The blue open circle shows the
estimated properties ($T_{\rm eff}\sim 4500$ K and luminosity $L\sim
3.25\times 10^{-3}\,L_{\odot}$) of the companion by Pan23. The solid black
line shows the $T_{\rm eff}$ and $L$ constraint from the $M_{\rm F606W}$
derived from the HST optical detection, and the black dashed lines correspond
to the 1$\sigma$ range.

##  3 Discussion & Summary

With the detection only in a single band (F606W), it is not possible to
independently derive COM-M71E’s physical properties from optical.
Nevertheless, we first make a comparison with the expected properties ($T_{\rm
eff}\sim 4500$ K and $L\sim 3.25\times 10^{-3}\,L_{\odot}$) from Pan23’s
preferred binary-evolution model. We use SYNPHOT function of MAAT (Ofek, 2014)
to estimate the bolometric correction in F606W for a blackbody at $T_{\rm
eff}=4500$ K, and find that the expected absolute magnitude is $M\sim 11.1$.
Given the potential theoretical uncertainties, it is broadly consistent with
$M_{\rm F606W}=11.7\pm 0.3$ from observation, supporting Pan23’s
interpretation that M71E is at an intermediate evolutionary stage between a RB
and a BW. Under the assumption of this model, we can estimate the physical
parameters of COM-M71E using our optical detection. By combining the
theoretical $\log(T_{\rm eff})-\log(L)$ track of the companion (blue solid
line in Figure 2) and the F606W photometric constraints (black lines n Figure
2), we obtain a blackbody radius $R_{\rm c,BB}\sim 0.08\,R_{\sun}$. Following
Pan23, we subsequently estimate its mass $M_{\rm c}\sim 0.09\,M_{\odot}$ based
on the mass-radius relations of Burrows et al. (1993), and thus $R_{\rm c,BB}$
is smaller than the Roche lobe radius $\sim 0.1\,R_{\odot}$. We note some
caveats regarding this interpretation. First of all, in Pan23’s simulation,
the BD donor’s surface is He-rich, whereas Burrows et al. (1993)’s BD models
assume solar abundance. Second, in this scenario, the companion was initially
H-rich. While its orbital period ($P\approx 53$ min) is larger than the
theoretically allowed minimum period of $P_{\rm min}\approx 37$ min for a
generic H-rich BD (Rappaport et al., 2021), it is smaller than the minimum
value $P_{\rm min}$ either predicted by binary evolution $(P_{\rm
min,theory}\sim 70$ min; Kolb & Baraffe 1999) or empirically seen in
cataclysmic variables $(P_{\rm min,observed}\sim 80$ min; Knigge 2006;
Gänsicke et al. 2009) for mass transfer from such a donor. Lastly, reaching a
mass of $M_{\rm c}\sim 0.09\,M_{\odot}$ requires an orbital inclination of
$\lesssim 5$ deg, which has a small geometric probability of $<1\%$ for
randomly distributed orbital orientations.

Pan23 disfavor WD companions because a plausible WD with the minimum mass of
$\sim 0.16\,M_{\odot}$ would need a finely-tuned face-on orbit with a very
small geometric probability ($<0.3\%$), and they also argue such a system is
challenging to form from theory. Such a minimum mass is based on the observed
low-mass WD population by Brown et al. (2016). A WD with $\approx
0.16\,M_{\odot}$ needs to have a temperature $\approx 10^{4}$K (typical for
WDs in Brown et al. 2016) to match the observed F606W flux. A minimum WD mass
of $\sim 0.16\,M_{\odot}$ is also consistent with the classic evolution
pathway by Tauris & Savonije (1999), but it corresponds to a minimum period
$\sim 1$ d, much larger than the observed $P\approx 53$ min. Alternatively, it
may be possible to reach such a short period from a partially evolved donor
star at the onset of mass transfer, and the orbital period shrinks during the
mass transfer, eventually leading to mass loss from the WD core (see, e.g.
Nelson et al., 1986; Deloye & Bildsten, 2003). Such a scenario may explain the
extremely low-mass He/C-rich degenerate companions found in some MSP systems
(e.g., Bailes et al. 2011; Romani et al. 2012). We briefly discuss the
possibility that M71E is an ultra-low-mass He-rich WD. We examine a fiducial
He WD with mass of $M_{\rm WD}=0.02\,M_{\odot}$ corresponding to a prior
likely inclination of $i\approx 30^{\circ}$, with a radius of $R_{\rm
WD}\approx 0.04\,R_{\odot}$ at $T_{\rm WD}\approx 6000$ K according to Deloye
& Bildsten (2003), it can match the observed F606W flux. Multi-band
photometric measurements would permit a temperature estimate and could thereby
test the WD scenario. Furthermore, spectroscopic observation can constrain its
chemical composition (H/He/C) and provide further clues to its formation
pathway.

Spider pulsar companions can be heated by PSR irradiation and can also be
significantly tidally distorted, inducing flux changes over its orbital
period. Known optical counterparts of spider pulsars can show significant
photometric variations modulated by the orbital period, usually with maxima at
inferior conjunctions (orbital phase $\phi=0.75$) and minima at superior
conjunctions ($\phi=0.25$). For example, the optical light curves of COM-M71A
(a BW companion in M71) show periodic variations with a full amplitude $\sim
3$ mag (Cadelano et al., 2015). Using Pan23’s pulsar timing parameters (Time
of Ascending Node at MJD $=58829.26006$ in solar system’s barycentric
reference frame and orbital period $P=0.0370398638$ d), our HST observation
(adding a -2.7 min correction to transform into the same reference frame as
pulsar timing) is at $\phi=0.67$, which is close to the expected maximum. The
light-curve amplitude is a function of orbital inclination (Pallanca et al.,
2014), and if the orbit is nearly face-on as expected by P23, the M71E system
would have much smaller variations than the eclipsing systems.

In conclusion, the HST optical observation of COM-M71A broadly agree with the
stripped dwarf model favored by Pan23. Future HST multi-band photometric
monitoring and JWST spectroscopic observations will be desirable for testing
various scenarios and making physical characterizations.

## Acknowledgment

We are grateful to the anonymous reviewer for insightful comments and
suggestions on the theoretical interpretation. We thank Jay Anderson, Hailiang
Chen, Andy Gould, Kejia Lee, Zhichen Pan and Bing Zhang for helpful
discussions. This work is supported by the science research grants from the
China Manned Space Project with No. CMS-CSST-2021-B12, the National Natural
Science Foundation of China (Grant No. 12133005) and the New Cornerstone
Science Foundation through the XPLORER PRIZE. This research is based on
observations made with the NASA/ESA Hubble Space Telescope obtained from the
Space Telescope Science Institute, which is operated by the Association of
Universities for Research in Astronomy, Inc., under NASA contract NAS 5-26555.
These observations are associated with program 16871. This work has made use
of data from the European Space Agency (ESA) mission Gaia
(https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and
Analysis Consortium (DPAC,
https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has
been provided by national institutions, in particular the institutions
participating in the Gaia Multilateral Agreement. This research has made use
of the VizieR catalogue access tool, CDS, Strasbourg, France.

## References

  * Alpar et al. (1982) Alpar, M. A., Cheng, A. F., Ruderman, M. A., & Shaham, J. 1982, Nature, 300, 728, doi: 10.1038/300728a0
  * Anderson (2022) Anderson, J. 2022, in American Astronomical Society Meeting Abstracts, Vol. 54, American Astronomical Society Meeting Abstracts, 206.01
  * Anderson et al. (2021) Anderson, J., Chiaberge, M., & Ryon, J. 2021, Astrometric Calibration of CTE, HST Proposal. Cycle 29, ID. #16871
  * Anderson & King (2000) Anderson, J., & King, I. R. 2000, PASP, 112, 1360, doi: 10.1086/316632
  * Bailes et al. (2011) Bailes, M., Bates, S. D., Bhalerao, V., et al. 2011, Science, 333, 1717, doi: 10.1126/science.1208890
  * Baumgardt & Vasiliev (2021) Baumgardt, H., & Vasiliev, E. 2021, MNRAS, 505, 5957, doi: 10.1093/mnras/stab1474
  * Benvenuto et al. (2014) Benvenuto, O. G., De Vito, M. A., & Horvath, J. E. 2014, ApJ, 786, L7, doi: 10.1088/2041-8205/786/1/L7
  * Bhattacharya & van den Heuvel (1991) Bhattacharya, D., & van den Heuvel, E. P. J. 1991, Phys. Rep., 203, 1, doi: 10.1016/0370-1573(91)90064-S
  * Brown et al. (2016) Brown, W. R., Gianninas, A., Kilic, M., Kenyon, S. J., & Allende Prieto, C. 2016, ApJ, 818, 155, doi: 10.3847/0004-637X/818/2/155
  * Burrows et al. (1993) Burrows, A., Hubbard, W. B., Saumon, D., & Lunine, J. I. 1993, ApJ, 406, 158, doi: 10.1086/172427
  * Cadelano et al. (2015) Cadelano, M., Pallanca, C., Ferraro, F. R., et al. 2015, ApJ, 807, 91, doi: 10.1088/0004-637X/807/1/91
  * Chen et al. (2013) Chen, H.-L., Chen, X., Tauris, T. M., & Han, Z. 2013, ApJ, 775, 27, doi: 10.1088/0004-637X/775/1/27
  * Deloye & Bildsten (2003) Deloye, C. J., & Bildsten, L. 2003, ApJ, 598, 1217, doi: 10.1086/379063
  * Dotter et al. (2010) Dotter, A., Sarajedini, A., Anderson, J., et al. 2010, ApJ, 708, 698, doi: 10.1088/0004-637X/708/1/698
  * Eggleton (1983) Eggleton, P. P. 1983, ApJ, 268, 368, doi: 10.1086/160960
  * Gaia Collaboration et al. (2016) Gaia Collaboration, Prusti, T., de Bruijne, J. H. J., et al. 2016, A&A, 595, A1, doi: 10.1051/0004-6361/201629272
  * Gaia Collaboration et al. (2023) Gaia Collaboration, Vallenari, A., Brown, A. G. A., et al. 2023, A&A, 674, A1, doi: 10.1051/0004-6361/202243940
  * Gänsicke et al. (2009) Gänsicke, B. T., Dillon, M., Southworth, J., et al. 2009, MNRAS, 397, 2170, doi: 10.1111/j.1365-2966.2009.15126.x
  * Ginzburg & Quataert (2020) Ginzburg, S., & Quataert, E. 2020, MNRAS, 495, 3656, doi: 10.1093/mnras/staa1304
  * Han et al. (2021) Han, J. L., Wang, C., Wang, P. F., et al. 2021, Research in Astronomy and Astrophysics, 21, 107, doi: 10.1088/1674-4527/21/5/107
  * Jia & Li (2015) Jia, K., & Li, X.-D. 2015, ApJ, 814, 74, doi: 10.1088/0004-637X/814/1/74
  * King et al. (2005) King, A. R., Beer, M. E., Rolfe, D. J., Schenker, K., & Skipp, J. M. 2005, MNRAS, 358, 1501, doi: 10.1111/j.1365-2966.2005.08877.x
  * King et al. (2003) King, A. R., Davies, M. B., & Beer, M. E. 2003, MNRAS, 345, 678, doi: 10.1046/j.1365-8711.2003.06990.x
  * Knigge (2006) Knigge, C. 2006, MNRAS, 373, 484, doi: 10.1111/j.1365-2966.2006.11096.x
  * Kolb & Baraffe (1999) Kolb, U., & Baraffe, I. 1999, MNRAS, 309, 1034, doi: 10.1046/j.1365-8711.1999.02926.x
  * Nelson et al. (1986) Nelson, L. A., Rappaport, S. A., & Joss, P. C. 1986, ApJ, 311, 226, doi: 10.1086/164767
  * Ofek (2014) Ofek, E. O. 2014, MAAT: MATLAB Astronomy and Astrophysics Toolbox, Astrophysics Source Code Library, record ascl:1407.005.  http://ascl.net/1407.005
  * Pallanca et al. (2014) Pallanca, C., Ransom, S. M., Ferraro, F. R., et al. 2014, ApJ, 795, 29, doi: 10.1088/0004-637X/795/1/29
  * Pan et al. (2023) Pan, Z., Lu, J. G., Jiang, P., et al. 2023, Nature, doi: 10.1038/s41586-023-06308-w
  * Rappaport et al. (2021) Rappaport, S., Vanderburg, A., Schwab, J., & Nelson, L. 2021, ApJ, 913, 118, doi: 10.3847/1538-4357/abf7b0
  * Roberts (2013) Roberts, M. S. E. 2013, in Neutron Stars and Pulsars: Challenges and Opportunities after 80 years, ed. J. van Leeuwen, Vol. 291, 127–132, doi: 10.1017/S174392131202337X
  * Romani et al. (2012) Romani, R. W., Filippenko, A. V., Silverman, J. M., et al. 2012, ApJ, 760, L36, doi: 10.1088/2041-8205/760/2/L36
  * Sarajedini et al. (2007) Sarajedini, A., Bedin, L. R., Chaboyer, B., et al. 2007, AJ, 133, 1658, doi: 10.1086/511979
  * Sirianni et al. (2005) Sirianni, M., Jee, M. J., Benítez, N., et al. 2005, PASP, 117, 1049, doi: 10.1086/444553
  * Tauris & Savonije (1999) Tauris, T. M., & Savonije, G. J. 1999, A&A, 350, 928, doi: 10.48550/arXiv.astro-ph/9909147