# SURDS: Self-Supervised Attention-guided Reconstruction and Dual Triplet Loss
for Writer Independent Offline Signature Verification

Soumitri Chattopadhyay
Saumik Bhattacharya  Dept. of Information Technology
Jadavpur University, Kolkata
Email<EMAIL_ADDRESS>Dept. of Electronics and Electrical
Communication Engineering
Indian Institute of Technology Kharagpur
Email<EMAIL_ADDRESS>Siladittya Manna
Umapada Pal  Computer Vision and Pattern Recognition Unit
Indian Statistical Institute, Kolkata
Email<EMAIL_ADDRESS>Computer Vision and Pattern Recognition Unit
Indian Statistical Institute, Kolkata
Email<EMAIL_ADDRESS>
###### Abstract

Offline Signature Verification (OSV) is a fundamental biometric task across
various forensic, commercial and legal applications. The underlying task at
hand is to carefully model fine-grained features of the signatures to
distinguish between genuine and forged ones, which differ only in minute
deformities. This makes OSV more challenging compared to other verification
problems. In this work, we propose a two-stage deep learning framework that
leverages self-supervised representation learning as well as metric learning
for writer-independent OSV. First, we train an image reconstruction network
using an encoder-decoder architecture that is augmented by a 2D spatial
attention mechanism using signature image patches. Next, the trained encoder
backbone is fine-tuned with a projector head using a supervised metric
learning framework, whose objective is to optimize a novel dual triplet loss
by sampling negative samples from both within the same writer class as well as
from other writers. The intuition behind this is to ensure that a signature
sample lies closer to its positive counterpart compared to negative samples
from both intra-writer and cross-writer sets. This results in robust
discriminative learning of the embedding space. To the best of our knowledge,
this is the first work of using self-supervised learning frameworks for OSV.
The proposed two-stage framework has been evaluated on two publicly available
offline signature datasets and compared with various state-of-the-art methods.
It is noted that the proposed method provided promising results outperforming
several existing pieces of work.

##  I Introduction

Handwritten signature is the most common and popular mode of biometric
identification that has been used over decades across various aspects of human
activity. From monetary transaction cheques to legal documents to offline
application forms, signature has been the de facto standard method for
document validation and is used almost at a daily basis. Thus, verification of
a signature is a very crucial task so as to accurately assert the original
owner of a signature sample, since fraudulence is widespread and often, highly
skilled. Such verification is heavily time consuming as well as prone to human
error, thereby required to be automated by pattern recognition methods.

Figure 1: Schematic of the intuition behind using the dual triplet objective
to fine-tune the pre-trained self-supervision model. More details have been
provided in Section III-C.

As such, signatures maybe categorised as (1) online and (2) offline, depending
on the modality it is obtained. Online signatures contain timestep-wise
coordinate information of a signature, usually obtained from an electronic
writing pad. Since such devices may be associated with sensors that can record
several parameters such as pen pressure, time taken per character and so on,
which are unique to a particular writer, verifying online signatures is a
fairly easy task. On the other hand, offline signatures are 2D images that are
captured electronically or manually before being processed for verification.
Such images contain no explicit information that can define a writing
behaviour, thereby requiring minute modeling of fine-grained features such as
character strokes and handwriting style. Thus, purely vision-based signature
verification is a challenging task and thus, has been widely studied over the
last few decades.

Broadly speaking, two approaches have been devised for OSV, they being (1)
writer dependent and (2) writer independent. In case of writer dependent OSV,
a classifier is trained for every writer class, which poses scalability
problems as they require the system to be re-trained each time a new user is
to be registered, as well as requires a large number of signature samples from
every user which is not easily available at one go. However, writer
independent OSV systems are more focused towards finding intrinsic differences
between genuine and forged samples such that they can be easily extended to
new users off-the-shelf without the need of re-training the system every time.
Thus, it is intuitive that the latter scenario is more close to the real-world
that can be easily up-scaled to industry standards.

In this research, we aim to develop a novel framework that leverages self-
supervised representation learning followed by supervised fine-tuning for
writer independent OSV. Although self-supervision has been employed to various
domains of computer vision, their application to offline handwritten data has
been limited. This might be due to the fact that existing self-supervised
approaches such as image colorization [1], super-resolution [2] or even
contrastive learning-based methods [3, 4] are designed for pixel perfect
natural or biomedical images [5] and thus, may not be suitable for sparse
binarized data [6]. In fact, to the best of our knowledge there has been no
prior work on offline signature verification that has employed self-
supervision strategy, which has been a driving motivation for our research.

In particular, we propose an encoder-decoder network for self-supervised image
reconstruction that is augmented by a patch-wise 2D attention module. The
attention mechanism allows the model to focus on visually important regions in
the entire image in a patch-wise manner. The intuition behind this mechanism
is to identify fine-grained local information from the signature samples.
Next, we incorporate a metric learning approach for fine-tuning the pre-
trained encoder by formulating an objective that separates negative samples
from both within the same writer as well as samples from different writers.
Specifically, we use a dual triplet loss module that is computed using
representations from both intra-writer and cross-writer negative samples with
respect to an anchor-positive pair of intra-writer signatures. The cross-
writer samples are implicitly negative samples as they come from different
writers. The intuition has been schematically depicted in Figure 1. This
allows a superior separation of genuine and forgeries, which has been
empirically shown while comparing with methods that use a single metric
learning objective [7]. Finally, the frozen encoder network is taken off-the-
shelf and used as a feature extractor for evaluation purpose. We have
evaluated the proposed model on two publicly available signature datasets [8],
where it has yielded results that are better or at par with other state-of-
the-art works in literature.

To summarize, the main contributions of the presented work are as follows:

  * •

A novel two-stage framework, named SURDS, has been proposed for writer
independent OSV.

  * •

A self-supervised attention-guided reconstruction-based pre-training strategy
has been proposed for learning fine-grained representations from signature
images.

  * •

A dual triplet loss-based fine-tuning is also proposed to tackle both intra-
writer and cross-writer negative samples.

  * •

To the best of our knowledge, this is the first work to employ self supervised
learning concept in the context of OSV.

  * •

The proposed framework outperforms adapted baseline competitors and also
surpasses several state-of-the-art works on two publicly available offine
signature datasets.

The rest of the paper is organized as follows: Section II discusses some of
the works in literature pertaining to OSV and self-supervised learning
approaches; Section III presents a detailed explanation of the proposed
framework for signature verification from offline images; Section IV describes
the implementation procedure and experimental protocol used in this work;
Section V analyses the empirical performance of the proposed model against
adapted baselines and state-of-the-art works; and finally, Section VI
concludes the findings of the proposed research.

##  II Related Works

###  II-A Offline Signature Verification

OSV has been well studied in recent decades, with the surveys [9, 10]
providing comprehensive review of the literature. Although classically
approached using handcrafted feature extraction and analysis [8, 11, 12], deep
learning-based methods have also been introduced in recent years, such as the
works presented in [7] that introduced contrastive loss-based Siamese CNNs for
signature verification; [13] which proposed a hierarchical CNN for learning
features from signature patches; and [14] that employed a hybrid two-channel
CNN for OSV. More recent works include region-based deep metric learning [15],
neuromotor inspired framework [16], Siamese networks followed by support
vector classifier [17], capsule networks-based model [18] and a graph neural
network approach [19]. Typically, OSV has been approached using fully
supervised methods only. In this work, we try to develop a self-supervised
pre-training method for writer independent signature verification.

###  II-B Self-Supervised Learning

Self-supervised learning [20, 3] has emerged as a powerful paradigm that seeks
to learn meaningful representations from unlabelled data so as to recognise
underlying patterns without the need for explicit annotation. This largely
alleviates the problem of scarcity of labelled data for specialised tasks,
since the unsupervised pre-trained representations are transferable to the
downstream task at hand. Classical means of self-supervised learning include
generative modelling [21, 22] and contrastive learning approaches [4, 23].
Particularly in the context of computer vision, self-supervision has been
explored by several pretext tasks for representation learning such as image
colorization [1], inpainting [24], super-resolution [2] and solving jigsaw
puzzles [25], among others. However, such methods are designed keeping pixel
perfect natural images in mind and hence are not applicable to sparse
greyscale signatures. One of the few works designed for handwritten data
employed cross-modal translation [6] between vector and raster space, which
however is unsuitable for our purpose since it requires the availability of
paired offline images and online coordinates. To this end, we propose a first-
of-its-kind self-supervised pre-training strategy for OSV that uses a patch-
wise attention guided reconstruction network for learning visual
representations from signature images.

##  III Methodology

###  III-A Preprocessing of Signature Images

Since raw signature images contain a lot of background pixels, directly
dividing them into patches would lead to a lot of redundant information and
thereby lead to poor reconstruction. To alleviate this, we devise a
preprocessing technique to ensure that the signatures are cropped to a tight
bound so as to minimize background pixels. We do so by first computing the
center of mass of the raw binarized image, after which we traverse through
vertical and horizontal axes on either side of the center of mass coordinates
to find the location of the last foreground pixel, at which we crop the image.
Once cropped to a tight bound, the images are resized to $256\times 256$
pixels and normalized before being fed into the encoder network. Further, this
resized image is divided into non-overlapping patches of $64\times 64$ pixels,
thereby obtaining $16$ patches for each image. A figure illustrating the
process has been shown in Figure 2.

Figure 2: Illustration of the preprocessing steps applied on the raw signature
images before feeding them into the encoder network.

###  III-B Self-Supervised Pre-training

Figure 3: Graphical representation of the proposed self-supervised
reconstruction-based pre-training pipeline for offline signature verification.

The goal of any self-supervised pipeline [1, 24] is to design a pre-training
strategy by which a model learns meaningful representations that can be used
off-the-shelf or fine-tuned for a downstream task. For signature images, we
set the goal to learn representations that contain local region level
information. To achieve this goal, we have designed a patch-wise attention-
guided image reconstruction network. Thus, it is evident that the core of the
pipeline would comprise an encoder-decoder architecture, the details of which
have been put forth in the following section.

####  III-B1 Convolutional Encoder-Decoder Architecture

For the convolutional encoder, we used the fairly straightforward ResNet-18
[26] model, which takes in a signature image $I\in\mathcal{R}^{H\times W\times
3}$ and outputs a convolutional feature map
$\mathcal{C}\in\mathcal{R}^{h\times w\times 512}$ having $512$ channels. While
feeding the full image into the network we perform a global average pooling on
the obtained feature map to get a feature vector of dimension $512$, whereas
the pooling operation is absent when the signature patches are fed. This is
due to the fact that we intend to perform a 2D attention between the patch
feature maps and the global feature vector, the details of which have been put
in the subsequent section.

For the decoder network, we have used the state-of-the-art U-Net [27] decoder
architecture. The decoder takes the attention-enriched $512$-dimensional
feature vector as input and performs a series of upsampling operations
comprising transpose convolutions to finally output a tensor
$Y\in\mathcal{R}^{H\times W\times 3}$ i.e. the dimensions of the original
signature image $I$. We have excluded the skip connections from the encoder
which are present in the original paper [27], where the authors aimed to
perform semantic segmentation. We formulate the reconstruction loss as the
simple mean squared error between the original and reconstructed images, which
thereby constitutes the self-supervision task objective given by Equation 1.

$\mathcal{L}_{SSL}=\frac{1}{H\times W}\sum_{i=1}^{H}\sum_{j=1}^{W}\lvert
I_{ij}-Y_{ij}\rvert^{2}$ (1)

####  III-B2 Patch-wise 2D Attention Mechanism

Figure 4: Schematic diagram showing the 2D attention mechanism employed in the
proposed study. The “broadcast operation” refers to that of PyTorch [28]
framework by which tensors of smaller dimensions are vectorized for
mathematical operations with a bigger dimensional tensor.

This is the key module in the self-supervised training pipeline that enforces
focusing on visually important parts of the signature in a region-wise
fashion. We use the convolution feature maps $\mathcal{C}$ obtained from the
signature patches and the feature vector $\mathcal{Z}$ obtained from the image
$I$ using the same weight shared encoder. The module outputs an attention
enriched context vector $\mathcal{H}$ using the equations given below:

$\begin{split}&\mathcal{M}=tanh(\mathcal{W}_{f}\circledast\mathcal{C}+\mathcal{W}_{h}\mathcal{Z})\\\
&\mathcal{A}_{ij}=softmax(\mathcal{W}_{att}\mathcal{M}_{ij})\\\
&{\color[rgb]{0,0,0}\mathcal{H}_{ij}=\sum_{k}\mathcal{A}^{T}_{ik}\cdot\mathcal{C}_{kj}}\end{split}$
(2)

Here, $\mathcal{W}_{f}$, $\mathcal{W}_{h}$ and $\mathcal{W}_{att}$ are
learnable weights, while the “$\circledast$” denotes the convolution
operation. The attention weight at $(i,j)$ is computed by $\mathcal{A}_{ij}$,
while $\mathcal{H}$ denotes the attention enriched context feature vector
obtained. For every patch feature map $\mathcal{C}_{k}$ we obtain a context
vector $\mathcal{H}_{k}$, which are then accumulated and a global average
pooling operation is performed to obtain a single feature vector of $512$
dimensions. This vector is then passed into the decoder for reconstruction.

###  III-C Metric Learning using Dual Triplet Loss

Once the pre-training using self-supervision is completed, we take the trained
ResNet-18 encoder and attach a projector head comprising two embedding layers
of $512$ dimensions each. This network is then fine-tuned using a metric
learning objective that distances a given query signature from both intra-
writer and cross-writer negatives, while bringing a positive signature closer
to itself (please refer to Figure 1 for a graphical form of the described
intuition). To do so, we have used the triplet loss function, a widely used
loss objective in deep metric learning that takes the representation of a
query (anchor, $A$) sample and minimizes its distance with a positive class
sample ($P$) and simultaneously maximizes its distance with a negative class
sample ($N$). The individual losses have been defined in Equation 3.

$\begin{split}&\mathcal{L}_{IW}=\max\\{0,\delta(z_{A},z_{P})-\delta(z_{A},z_{N_{IW}})+\mu\\}\\\
&\mathcal{L}_{CW}=\max\\{0,\delta(z_{A},z_{P})-\delta(z_{A},z_{N_{CW}})+\mu\\}\end{split}$
(3)

Where, $\mu>0$ is a hyperparameter used to avoid trivial solution, while
$\delta(z_{i},z_{j})=\lVert z_{i}-z_{j}\rVert_{2}$.

Combining the intra-writer and cross-writer losses, we define our dual triplet
objective as:

$\mathcal{L}_{DT}=\mathcal{L}_{IW}+\mathcal{L}_{CW}$ (4)

###  III-D Downstream Evaluation

After the training phase is completed, we take the frozen backbone comprising
the ResNet encoder and projector and use it off-the-shelf for feature
extraction from unseen signature images for downstream evaluation. We have
used the evaluation protocol described by Dey et al. [7] for writer
independent signature verification. Specifically, we have used a threshold
$\tau$ that determines if a pair of signature $(i,j)$ is genuine or forged
depending on a distance metric $\mathcal{D}(z_{i},z_{j})$ between their
embeddings, as given by Equation 5. However, genuineness/falsification is a
relative decision i.e. it is decided using a ground truth reference. For this
purpose, we randomly select out $8$ genuine samples from each of the unseen
writers and consider them as reference signatures, and the mean of the
distance metric between each of the references with the queried sample is
considered for comparison with the threshold. The remaining samples are then
used for evaluation.

$\mathcal{D}(z_{i},z_{j})=\lVert z_{i}-z_{j}\rVert$ (5)

##  IV Experiments

###  IV-A Datasets

We have used two publicly available offline signature datasets, BHSig260
Bengali and Hindi [8] for evaluating our proposed two-stage pipeline. The
Bengali signature dataset comprises signatures of 100 writers, while the Hindi
dataset consists of signatures coming from 160 writers. For each writer, there
are with 24 genuine signature samples and 30 skilled forgeries. We have split
the datasets randomly in 0.7/0.3 ratio for training/testing i.e. 70/30 writers
for the Bengali dataset and 112/48 writers for the Hindi dataset respectively.
Since our premise is writer independent signature verification, we have
ensured that the writers _do not_ overlap between the respective sets.

###  IV-B Implementation

We have implemented our proposed framework using the PyTorch framework [28] on
a 8GB Nvidia GeForce RTX 2080 GPU. For the CNN encoder used in self-supervised
pre-training, we have used a randomly initialized ResNet-18 [26] network
excluding the linear layers, which is then re-used off-the-shelf for the
metric learning phase after pre-training is completed. The optimizer used for
the pre-training task is SGD [29] with an initial learning rate of $0.1$. For
the fine-tuning phase, we have used the RMSprop optimizer with initial
learning rates for the encoder and the projector being $0.005$ and $1.0$
respectively and weight decay factor of $0.0005$. For both the phases, the
respective models have been trained for 200 epochs, and a Linear Warmup Cosine
Annealing scheduler has been used for learning rate decay. Prior to being fed
into the encoder, the signature images have been preprocessed using the steps
described in Section III-A.

###  IV-C Evaluation Metrics

For signature pairs $(i,j)$ let $\mathcal{G}$ and $\mathcal{G^{\prime}}$
denote sets of genuine-genuine and genuine-forged signature pairs
respectively. Then, at a given threshold $\tau$ we can define the set of true
positives ($T_{P}$), true negatives ($T_{N}$), false positives ($F_{P}$) and
false negatives ($F_{N}$) as:

$\begin{split}&T_{P}(\tau)=\\{(i,j)\in\mathcal{G},\mathcal{D}(z_{i},z_{j})\leq\tau\\}\\\
&T_{N}(\tau)=\\{(i,j)\in\mathcal{G^{\prime}},\mathcal{D}(z_{i},z_{j})>\tau\\}\\\
&F_{P}(\tau)=\\{(i,j)\in\mathcal{G^{\prime}},\mathcal{D}(z_{i},z_{j})\leq\tau\\}\\\
&F_{N}(\tau)=\\{(i,j)\in\mathcal{G},\mathcal{D}(z_{i},z_{j})>\tau\\}\end{split}$
(6)

Accordingly, the respective rates may be defined as the ratio of set
cardinality to its superset cardinality, denoted by $n(\cdot)$.

$\begin{split}&TPR(\tau)=\frac{n(T_{P}(\tau))}{n(\mathcal{G})};\;\;\;\;TNR(\tau)=\frac{n(T_{N}(\tau))}{n(\mathcal{G^{\prime}})}\\\
&FPR(\tau)=\frac{n(F_{P}(\tau))}{n(\mathcal{G^{\prime}})};\;\;\;\;FNR(\tau)=\frac{n(F_{N}(\tau))}{n(\mathcal{G})}\end{split}$
(7)

Based on these, we define the three evaluation metrics used in this work:

  * •

Accuracy:  The accuracy is defined as the arithmetic mean of the true positive
and true negative rates. We report the maximum accuracy achieved by varying
the threshold $\tau$ from the minimum distance to maximum distance value at
intervals of $5e-5$.

$Accuracy=\max\limits_{\tau\in\mathcal{D}}\left\\{\frac{TPR(\tau)+TNR(\tau)}{2}\right\\}$
(8)

  * •

False Acceptance Rate (FAR):  FAR is defined as the ratio of forged samples
incorrectly determined to be genuine and thereby accepted (i.e. $FPR$),

  * •

False Rejection Rate (FRR):  FRR is defined as the ratio of genuine samples
incorrectly predicted to be forgeries and thereby rejected (i.e. $FNR$).

##  V Results and Analysis

###  V-A Ablation Study

We first compare our proposed framework with some ablated baseline competitor
models derived from our architecture:

  * •

RN-DTL:  We ablate the self-supervised pre-training process and train a
randomly initialized ResNet-18 encoder using the dual triplet loss based
supervised metric learning proposed in our study.

  * •

AE-DTL:  We ablate the attention mechanism and simply train a convolutional
autoencoder using the encoder-decoder architecture described in Section
III-B1, followed by fine-tuning of the encoder using our proposed dual triplet
loss module.

The results of the comparative study have been tabulated in Table I. From the
results, it is evident that self-supervised pre-training significantly boosts
the downstream signature verification performance.

TABLE I: Results obtained from the ablation study conducted on our proposed
framework.

Baseline |  BHSig260 Bengali [8] |  BHSig260 Hindi [8]
---|---|---
Accuracy (%) | FAR (%) | FRR (%) | Accuracy (%) | FAR (%) | FRR (%)
RN-DTL | 82.18 | 21.66 | 13.95 | 86.72 | 14.67 | 11.88
AE-DTL | 83.57 | 29.11 | 3.75 | 89.14 | 18.47 | 3.25
Proposed Method | 87.34 | 19.89 | 5.42 | 89.50 | 12.01 | 8.98

###  V-B Comparison with state-of-the-art

To the best of our knowledge, there is no prior work that uses self-supervised
pre-training for writer independent OSV, and hence, we compare our proposed
framework with some of the existing state-of-the-art supervised methods in
literature. The empirical results are shown in Table II. It is evident that
our proposed pipeline shows promising performance on both the datasets,
outperforming several of the existing works by large margins. Furthermore, the
FRR values obtained by the proposed method are significantly better than all
of the existing works. The empirical comparisons highlight the usefulness and
reliability of the proposed method for signature verification.

TABLE II: Comparison of our proposed framework with some state-of-the-art
works in literature.

Method |  BHSig260 Bengali [8] |  BHSig260 Hindi [8]
---|---|---
Accuracy (%) | FAR (%) | FRR (%) | Accuracy (%) | FAR (%) | FRR (%)
Pal et al. [8] | 66.18 | 33.82 | 33.82 | 75.53 | 24.47 | 24.47
Dey et al. [7] | 86.11 | 13.89 | 13.89 | 84.64 | 15.36 | 15.36
Dutta et al. [11] | 84.90 | 15.78 | 14.43 | 85.90 | 13.10 | 15.09
Alaei et al. [12] | – | 16.18 | 30.12 | – | 16.18 | 30.12
Rateria et al. [17] | 75.06 | 27.81 | 21.74 | 89.33 | 10.93 | 10.02
Proposed Method | 87.34 | 19.89 | 5.42 | 89.50 | 12.01 | 8.98

###  V-C Cross-Dataset Evaluation

Further, we extend our evaluations to a cross-dataset scenario i.e. training
on one dataset and evaluating on a different dataset. The aim of this
investigation is to find out the transferability of representations obtained
from signatures of one language to another. To do so, We performed the
experiments on two possible setups:

  * •

Setup-1:  We train only the self-supervision pipeline on one dataset and
perform the fine-tuning phase as well as evaluate on the other. This is much
like a fine-tuning based transfer learning setup. Since the overall setup
requires training samples from both datasets, we have maintained the
train/test split as described in Section IV-A. The rest of the framework
remains unchanged.

  * •

Setup-2:  We perform the entire training process, self-supervised training
followed by dual triplet loss-based fine-tuning, on one dataset. We then use
this fine-tuned model for inference on the downstream dataset (which contains
signatures from a different language). For this setup, we evaluate and report
our results on the entire dataset to be used for testing. Although we could
have as well trained our model on an entire dataset for this setup, we
nevertheless chose to follow the training set configuration as given in
Section IV-A so that our training process remains consistent with those
reported in Tables I and II.

The results are shown in Table III. We observe that in Setup-1 (i.e. fine-
tuning based transfer learning setup), the results obtained are almost at par
with state-of-the-art works as well as with the proposed method in an intra-
dataset setup (tabulated in Table II). For Setup-2, the performance
detoriation obtained is expected, since all datasets have some intrinsically
unique features which can be learnt only from intra-dataset training. Further,
it must also be noted that in Setup-2 we have not used the entire training
dataset and instead used only a fraction of it to maintain consistency with
our original experiments. Nevertheless, the cross-dataset transfer learning
results for Setup-2 are empirically better than those obtained by Dey et al.
[7], where the authors reported accuracies of 64.57% (train: Bengali, test:
Hindi) and 60.65% (train: Hindi, test: Bengali) on cross-dataset evaluation,
which highlights the effectiveness of our proposed signature representation
learning process that can be ably transferred across datasets.

TABLE III: Cross-dataset evaluation results obtained by the proposed method.

Train/Test | Setup-1 | Setup-2
---|---|---
Accuracy (%) | FAR (%) | FRR (%) | Accuracy (%) | FAR (%) | FRR (%)
Bengali/Hindi | 88.19 | 20.62 | 2.99 | 67.01 | 51.14 | 14.82
Hindi/Bengali | 86.14 | 17.11 | 10.01 | 72.26 | 46.46 | 9.00

###  V-D Cross-Script Evaluation

Finally, we also probe into a cross-script OSV setup, where the model trained
on a signature corpus of a given script is evaluated on a dataset of a
different script altogether. For this setup, we first train our proposed
framework separately on the Bengali and Hindi datasets (i.e. Indic script) [8]
used in this study and then test the trained models on the ICDAR-2011 Dutch
signature dataset [30]. The Dutch dataset has a pre-defined evaluation set
comprising query (1287) as well as reference (646) signatures coming from 54
writers. Here also, we have adhered to the training set distribution described
in Section IV-A. The results obtained are:

  * •

Train on Bengali:  Accuracy = $68.60\%$, FAR = $43.51\%$, FRR = $19.29\%$

  * •

Train on Hindi:  Accuracy = $76.65\%$, FAR = $40.06\%$, FRR = $6.63\%$

It is worth noting that our setup is mere transfer learning without any
explicit fine-tuning and thus, the model is completely unaware of any script
specific knowledge. Yet, our model shows promise in the evaluation, especially
with its low FRR values, which is commendable. The results show that our study
has the potential to be explored in developing a script-independent universal
signature verification framework.

##  VI Conclusion and Future Work

In this study, we have proposed a novel two-stage pipeline for offline
signature verification that employs a self-supervised reconstruction network
for pre-training, which is then fine-tuned using a metric learning objective.
The reconstruction model is ably augmented using a patch-wise 2D attention
module which enforces the network to focus on local signature patches so as to
learn meaningful representations that can be transferred to the downstream
verification task. The fine-tuning phase uses a dual triplet objective that
effectively discriminates between genuine and forged signature
representations. We have shown by empirical results that our model shows
promising performance on publicly available benchmark signature datasets,
outperforming quite a few existing state-of-the-art works in literature.
Further, comparisons with baseline competitors highlight the usefulness of the
respective components of the proposed framework, and testing in challenging
cross-dataset and cross-script setups emphasize its application in a diverse
environment. Being first-of-its-kind, this work potentially paves the way for
developing self-supervised pre-training tasks for the domain of signature
verification, which is very important due to the bottleneck in availability of
annotated data in the real world. Furthermore, the proposed pre-training
pipeline may also be employed to other domains of computer vision, especially
where region-wise focusing is crucial. We intend to work further in these
directions in our future works.

## References

  * [1] R. Zhang, P. Isola, and A. A. Efros, “Colorful image colorization,” in _ECCV_ , 2016.
  * [2] C. Ledig, L. Theis, F. Huszár, J. Caballero, A. Cunningham, A. Acosta, A. Aitken, A. Tejani, J. Totz, Z. Wang _et al._ , “Photo-realistic single image super-resolution using a generative adversarial network,” in _CVPR_ , 2017.
  * [3] A. Jaiswal, A. R. Babu, M. Z. Zadeh, D. Banerjee, and F. Makedon, “A survey on contrastive self-supervised learning,” _Technologies_ , 2021.
  * [4] T. Chen, S. Kornblith, M. Norouzi, and G. Hinton, “A simple framework for contrastive learning of visual representations,” in _ICML_ , 2020.
  * [5] S. Shurrab and R. Duwairi, “Self-supervised learning methods and applications in medical imaging analysis: A survey,” _arXiv preprint arXiv:2109.08685_ , 2021.
  * [6] A. K. Bhunia, P. N. Chowdhury, Y. Yang, T. M. Hospedales, T. Xiang, and Y.-Z. Song, “Vectorization and rasterization: Self-supervised learning for sketch and handwriting,” in _CVPR_ , 2021.
  * [7] S. Dey, A. Dutta, J. I. Toledo, S. K. Ghosh, J. Lladós, and U. Pal, “Signet: Convolutional siamese network for writer independent offline signature verification,” _arXiv preprint arXiv:1707.02131_ , 2017.
  * [8] S. Pal, A. Alaei, U. Pal, and M. Blumenstein, “Performance of an off-line signature verification method based on texture features on a large indic-script signature dataset,” in _DAS_. IEEE, 2016.
  * [9] L. G. Hafemann, R. Sabourin, and L. S. Oliveira, “Offline handwritten signature verification—literature review,” in _IEEE IPTA_ , 2017.
  * [10] H. Kaur and M. Kumar, “Signature identification and verification techniques: state-of-the-art work,” _Journal of Ambient Intelligence and Humanized Computing_ , 2021.
  * [11] A. Dutta, U. Pal, and J. Lladós, “Compact correlated features for writer independent signature verification,” in _ICPR_. IEEE, 2016.
  * [12] A. Alaei, S. Pal, U. Pal, and M. Blumenstein, “An efficient signature verification method based on an interval symbolic representation and a fuzzy similarity measure,” _IEEE TIFS_ , 2017.
  * [13] S. Shariatmadari, S. Emadi, and Y. Akbari, “Patch-based offline signature verification using one-class hierarchical deep learning,” _IJDAR_ , 2019\.
  * [14] M. Berkay Yilmaz and K. Ozturk, “Hybrid user-independent and user-dependent offline signature verification with a two-channel cnn,” in _CVPRW_ , 2018\.
  * [15] L. Liu, L. Huang, F. Yin, and Y. Chen, “Offline signature verification using a region based deep metric learning network,” _Pattern Recognition_ , 2021\.
  * [16] M. Diaz, M. A. Ferrer, and R. Sabourin, “Approaching the intra-class variability in multi-script static signature evaluation,” in _ICPR_ , 2016\.
  * [17] A. Rateria and S. Agarwal, “Off-line signature verification through machine learning,” in _IEEE UPCON_ , 2018.
  * [18] E. Parcham, M. Ilbeygi, and M. Amini, “Cbcapsnet: A novel writer-independent offline signature verification model using a cnn-based architecture and capsule neural networks,” _Expert Systems with Applications_ , 2021.
  * [19] S. Roy, D. Sarkar, S. Malakar, and R. Sarkar, “Offline signature verification system: a graph neural network based approach,” _Journal of Ambient Intelligence and Humanized Computing_ , 2021.
  * [20] L. Jing and Y. Tian, “Self-supervised visual feature learning with deep neural networks: A survey,” _IEEE TPAMI_ , 2020.
  * [21] D. P. Kingma and M. Welling, “Auto-encoding variational bayes,” _arXiv preprint arXiv:1312.6114_ , 2013.
  * [22] I. Goodfellow, J. Pouget-Abadie, M. Mirza, B. Xu, D. Warde-Farley, S. Ozair, A. Courville, and Y. Bengio, “Generative adversarial nets,” in _NeurIPS_ , 2014.
  * [23] K. He, H. Fan, Y. Wu, S. Xie, and R. Girshick, “Momentum contrast for unsupervised visual representation learning,” in _CVPR_ , 2020.
  * [24] D. Pathak, P. Krahenbuhl, J. Donahue, T. Darrell, and A. A. Efros, “Context encoders: Feature learning by inpainting,” in _CVPR_ , 2016.
  * [25] M. Noroozi and P. Favaro, “Unsupervised learning of visual representations by solving jigsaw puzzles,” in _ECCV_ , 2016, pp. 69–84.
  * [26] K. He, X. Zhang, S. Ren, and J. Sun, “Deep residual learning for image recognition,” in _CVPR_ , 2016.
  * [27] O. Ronneberger, P. Fischer, and T. Brox, “U-net: Convolutional networks for biomedical image segmentation,” in _MICCAI_ , 2015.
  * [28] A. Paszke, S. Gross, F. Massa, A. Lerer, J. Bradbury, G. Chanan, T. Killeen, Z. Lin, N. Gimelshein, L. Antiga _et al._ , “Pytorch: An imperative style, high-performance deep learning library,” in _NeurIPS_ , 2019.
  * [29] I. Sutskever, J. Martens, G. Dahl, and G. Hinton, “On the importance of initialization and momentum in deep learning,” in _ICML_ , 2013.
  * [30] G. Alvarez, B. Sheffer, and M. Bryant, “Offline signature verification with convolutional neural networks,” _Technical report, Stanford University_ , 2016.


# Prospects for constraining twin stars with next-generation
gravitational-wave detectors

Philippe Landry  Canadian Institute for Theoretical Astrophysics, University
of Toronto, Toronto, Ontario M5S 3H8, Canada<EMAIL_ADDRESS>Kabir
Chakravarti  CEICO, FZU-Institute of Physics of the Czech Academy of Sciences,
Na Slovance 1999/2, 182 21 Prague 8, Czech Republic<EMAIL_ADDRESS>
###### Abstract

Neutron star equations of state with strong phase transitions may support twin
stars, hybrid and hadronic stars with the same mass but different tidal
deformabilities. The presence of twin stars in the population of merging
neutron stars produces distinctive gaps in the joint distribution of binary
tidal deformabilities and chirp masses. We analyze a simulated population of
binary neutron star mergers recovered with a network of next-generation (XG)
ground-based gravitational-wave detectors to determine how many observations
are needed to infer, or rule out, the existence of twin stars. Using a
hierarchical inference framework based on a simple parametric twin-star model,
we find that a single week of XG observations may suffice to detect a tidal
deformability difference of several hundred between twins and measure the mass
scale at which twins occur to within a few percent. For less pronounced twins,
XG observations will place a stringent upper bound on the tidal deformability
difference.

##  I Introduction

The phase structure of matter at the highest densities realized inside neutron
stars is an unsolved puzzle for nuclear physics. Above the nuclear saturation
density of $\rho_{\rm nuc}=2.8\times 10^{14}$ g/cm3, the nucleonic
constituents of ordinary matter are thought to give way to other fundamental
degrees of freedom, such as hyperons or deconfined quarks. If this phase
transition occurs at a density that prevails inside neutron stars, then the
heaviest of these compact objects may in fact be hybrid stars with exotic-
matter cores. The possible existence of a stable family of hybrid stars with
central densities greater than those of conventional hadronic neutron stars
has been envisioned and studied extensively from the theoretical point of view
Gerlach (1968); Kampfer (1981); Glendenning and Kettner (2000); Schertler  _et
al._ (2000); Alford  _et al._ (2013); Zdunik and Haensel (2013); Benić  _et
al._ (2015); Alford  _et al._ (2015); Han and Steiner (2019).

In many nuclear theory models, the transition to the high-density phase is of
first order, exhibiting a discontinuous first derivative of the baryon
chemical potential with respect to the baryon density $\rho$. The hadronic-
and exotic-matter phases may either remain separate and interface directly, or
they may be linked by an intermediate mixed phase. In the former scenario, the
so-called Maxwell construction, the baryon density experiences a discontinuity
$\Delta\rho$ at fixed pressure $p$. If $\Delta\rho$ is sufficiently large, the
dense matter equation of state may support stable twin stars: hadronic- and
hybrid-star counterparts with the same mass, but different compactness
Schaffner-Bielich  _et al._ (2002); Zacchi  _et al._ (2017); Alford and
Sedrakian (2017); Christian  _et al._ (2018).

Figure 1: Sequences of neutron star properties for selected twin-star-
supporting equations of state. We plot the mass-radius (top) and mass-tidal
deformability (bottom) relations. Dashed lines indicate unstable segments of
the sequence, where no neutron stars exist. Markers bookend the mass range
over which twins occur.

These twin-star configurations are visible in the mass-radius ($m$–$R$)
relations of Fig. 1, obtained by solving the Tolman-Oppenheimer-Volkoff (TOV)
equations Oppenheimer and Volkoff (1939); Tolman (1939) for selected equations
of state. In each relation, there is a small range of masses for which a
hybrid twin with a smaller radius and larger central density coexists
alongside a hadronic star. Similar morphology is visible in the mass-tidal
deformability ($m$–$\Lambda$) relations obtained by solving for a quadrupolar
tidal perturbation Flanagan and Hinderer (2008); Hinderer (2008); Landry and
Poisson (2014). The presence of radius or tidal deformability twins is a
smoking gun for the occurrence of a strong first-order phase transition in
neutron star matter.

Because twin stars provide such a clear signature in the $m$–$R$ (or
$m$–$\Lambda$) relation, their compatibility with existing observations of
neutron stars has been heavily scrutinized. References Paschalidis  _et al._
(2018); Montaña  _et al._ (2019); Christian  _et al._ (2019); Pang  _et al._
(2020); Wang  _et al._ (2022) demonstrated that parametric equations of state
supporting twin stars could be constructed to simultaneously satisfy
observational constraints on the tidal deformability from the compact binary
merger GW170817 Abbott  _et al._ (2017, 2019, 2018a) and bounds on the maximum
mass $M_{\rm TOV}$ from the heaviest pulsars discovered in radio surveys
Antoniadis  _et al._ (2013); Cromartie  _et al._ (2020); Fonseca  _et al._
(2021). Reference Christian and Schaffner-Bielich (2020) argued that the
radius measured for PSR J0030+0451 via X-ray pulse profile modeling Riley  _et
al._ (2019); Miller  _et al._ (2019) rules out twin stars associated with a
phase transition onset density below $1.7\rho_{\rm nuc}$. Twin stars were also
shown Christian and Schaffner-Bielich (2022) to be compatible with the
inferred radius of PSR J0740+6620 Riley  _et al._ (2021); Miller  _et al._
(2021) and Li  _et al._ (2021) with PREX-II’s measurement of the neutron skin
of 208Pb Adhikari  _et al._ (2021). Meanwhile, Ref. Christian and Schaffner-
Bielich (2021) suggested (contra Ref. Tsaloukidis  _et al._ (2022)) that a
maximum mass in excess of $2.2\,M_{\odot}$ would exclude twin stars across the
entire neutron star mass spectrum. Hence, although existing astronomical
observations substantially restrict the parameter space for twin stars, they
remain a live possibility.

Additional observations of neutron stars with the existing gravitational-wave
detector network of Advanced LIGO Aasi  _et al._ (2015), Virgo Acernese  _et
al._ (2015) and KAGRA Akutsu  _et al._ (2021) could, in principle, turn up a
serendipitous discovery of twin stars. However, this is unlikely given the
typically narrow mass range over which twins are supported and the high
resolution in $\Lambda$ required to identify them—not to mention the modest
expected binary neutron star detection rate Abbott  _et al._ (2018b); Colombo
_et al._ (2022); Patricelli  _et al._ (2022). Prospects for detecting twin
stars with next-generation (XG) ground-based gravitational-wave observatories
like Cosmic Explorer Evans  _et al._ (2021) and Einstein Telescope Maggiore
_et al._ (2020) are more promising, thanks to their ability to capture
virtually the complete merging neutron star population in the nearby Universe
Evans  _et al._ (2021); Borhanian and Sathyaprakash (2022). Although precisely
measuring the tidal deformability of individual neutron stars will remain a
challenging task,111Ref. Smith  _et al._ (2021) reports statistical
uncertainties of $O(100\%)$ in component tidal deformabilities $\Lambda_{1,2}$
even for a simulated XG binary neutron star merger at 40 Mpc with signal-to-
noise ratio 2400. twin-star-supporting equations of state give rise to a
distinctive distribution of binary tidal deformabilities $\tilde{\Lambda}$ vs
chirp masses $\mathcal{M}$ across the population Chatziioannou and Han (2020):
as shown in Fig. 2, when no twin stars are present, the distribution is
contiguous, while twins introduce gaps that make it disjoint. This population-
level signature can be used to infer the existence of twin stars even when
individual twin-star pairs are not identifiable. Here we show that XG
measurements of the binary tidal deformability distribution can detect the
existence of twin stars at 90% confidence (or with a log Bayes factor greater
than 6) with as few as $\sim 100$ observations.

To do so, we first construct a family of equations of state with first-order
phase transitions giving rise to stable twin stars, and simulate populations
of binary neutron star mergers recovered with different gravitational-wave
detector networks. We then collect the binary tidal deformability and chirp
mass measurements from each population. Heuristically, twin stars are detected
at the population level when the recovered distribution is disjoint, rather
than contiguous. To implement this test systematically, we develop a Bayesian
hierarchical inference framework with a parametric model for the $m$–$\Lambda$
relation, formulated in terms of a twin-star mass scale $M_{t}$ and the tidal
deformability difference $\Delta\Lambda$ between twins. This approach follows
Ref. Chatziioannou and Han (2020), which tackled the related problem of
inferring the radius difference between hybrid and hadronic stars with LIGO
and Virgo. We find that current gravitational-wave detectors, even at upgraded
“A+” sensitivities Abbott  _et al._ (2018b), cannot meaningfully constrain the
twin-star parameters. However, in the best-case scenario we consider, an XG
network can determine $M_{t}$ and $\Delta\Lambda$ to within 1% and 15%,
respectively, at 90% confidence after just one month of observations.

Figure 2: Joint binary tidal deformability and chirp mass distribution probed
by one month of XG observations for the twin-star-supporting SKI5_2009
equation of state, as compared with the purely hadronic SKI5 equation of
state. A uniform neutron star mass distribution is assumed. Note the gaps in
the distribution for SKI5_2009 indicated by the arrows, which signify the
existence of twin stars, compared to the contiguous SKI5 distribution. The
detected distribution is subject to statistical uncertainties and is scattered
by detector noise; for SKI5_2009, we show the recovered chirp masses and
binary tidal deformabilities in grey, with marker sizes inversely proportional
to the signal-to-noise ratio. Actual statistical uncertainties are $O(10^{2})$
in $\tilde{\Lambda}$ and $O(10^{-3})$ in $\mathcal{M}$.

##  II Twin-star equations of state

We construct equations of state giving rise to twin stars using the constant-
sound-speed formulation of Ref. Han and Steiner (2019). It consists of a low-
density hadronic equation of state, to which a constant-pressure first-order
phase transition segment of length (or “strength”) $\Delta\rho$ is connected
at onset density $\rho_{t}$. A high-density extension with constant sound
speed is appended to the end of the phase transition segment.

We demand that the equations of state satisfy existing constraints from
neutron star observations. If we suppose that the $\sim 1.4\,M_{\odot}$
neutron stars observed in GW170817 and as PSR J0030+0451 are purely hadronic,
then the low-density equation of state must satisfy the approximate
observational upper bounds $\Lambda_{1.4}\lesssim 580$ Abbott  _et al._
(2018a) and $R_{1.4}\lesssim 14$ km Miller  _et al._ (2019).222In fact,
because the analysis in Ref. Abbott  _et al._ (2018a) relies on equation-of-
state-insensitive relations that can’t accommodate twin stars, we allow for
some leeway in the first constraint. If, instead, we suppose that those
compact objects are hybrid stars, then the low-density portion of the equation
of state remains unconstrained. These considerations motivate the choice of
SKI272 ($R_{1.4}\approx 13.5$ km, $\Lambda_{1.4}\approx 650$) Agrawal  _et
al._ (2003) and SKI5 ($R_{1.4}\approx 14.5$ km, $\Lambda_{1.4}\approx 1000$)
Reinhard and Flocard (1995); Bertulani and Valencia (2019), respectively. Both
of these low-density equations of state have a symmetry energy slope
$L\lesssim 140$, in keeping with the implications of the PREX-II experiment
Reed  _et al._ (2021).

EOS | $\rho_{t}\;[\rho_{\rm nuc}]$ | $\Delta\rho\;[\rho_{t}]$ | $\Delta M_{t}\;[M_{\odot}]$ | $M_{\rm TOV}\;[M_{\odot}]$ |  $R_{1.4}\;$[km] | $\Lambda_{1.4}$ | $M_{t}\;[M_{\odot}]$ | $\Delta\Lambda$
---|---|---|---|---|---|---|---|---
SKI5_2006 | 2.0 | 0.6 | 0.01 | 2.28 | 13.5 | 520 | 1.36 | 340
SKI5_2009 | 2.0 | 0.9 | 0.07 | 2.08 | 12.2 | 240 | 1.32 | 940
SK272_2506 | 2.5 | 0.6 | 0.01 | 2.08 | 13.5 | 650 | 1.45 | 80
SK272_3505 | 3.5 | 0.5 | 0.03 | 1.94 | 13.5 | 650 | 1.91 | 40
Table 1: Twin-star-supporting equations of state used in this study. The
parameters $\rho_{t}$ and $\Delta\rho$ specifying the onset density and
strength of the first order phase transition are listed. The mass range
$\Delta M_{t}$ over which twins are supported, the maximum mass $M_{\rm TOV}$,
and the canonical radius $R_{1.4}$ and tidal deformability $\Lambda_{1.4}$ for
the corresponding sequence of neutron stars are also given. Additionally, we
list the twin-star mass scale $M_{t}$ and tidal deformability difference
$\Delta\Lambda$ that we seek to recover via hierarchical inference. The
equation of state’s name indicates the low-density hadronic model upon which
it is based.

For the high-density extension, we select a causal equation of state to
increase the odds that stable twin stars will be supported. To interpolate
between the high- and low-density regimes, we explore different first-order
phase transition segments by selecting various combinations of onset density
$\rho_{t}\in[1,4]\rho_{\rm nuc}$ and strength $\Delta\rho\in(0,2]\rho_{t}$.
For use in our study, we retain a representative subset of those equations of
state that give rise to stable twin stars while simultaneously satisfying the
approximate mass and radius constraints from electromagnetic observations of
PSR J0740+6620: $M_{\rm TOV}\gtrsim 2.0\,M_{\odot}$, $R_{2.0}\gtrsim 12$ km
Miller  _et al._ (2021). The selected equations of state, their phase
transition parameters and their associated neutron star observables are listed
in Table 1. Their $m$–$R$ and $m$–$\Lambda$ relations are plotted in Fig. 1.
Of particular relevance for our analysis is the presence of an unstable
segment in the $m$–$\Lambda$ relation that connects the hadronic and hybrid
branches at a mass scale corresponding to $\rho_{t}$.

##  III Binary neutron star population and tidal deformability distribution

To determine when twin stars are identifiable in the binary neutron star
population, we generate simulated distributions of masses and binary tidal
deformabilities. We prescribe a neutron star mass model, select an equation of
state, sample a realization of the astrophysical population of binary neutron
star mergers, and determine which events are detected by computing their
optimal signal-to-noise ratio with respect to the detector networks of
interest.

Inspired by studies of the gravitational-wave population to date Landry and
Read (2021); Abbott  _et al._ (2021), our fiducial neutron star mass
distribution is uniform for $m\in[1\,M_{\odot},M_{\rm TOV}]$. We assume that
both components of a neutron star binary are drawn from this common mass
distribution, and that they pair randomly. We distribute the sources
isotropically on the sky and according to a Madau-Dickinson star formation
rate Madau and Dickinson (2014) in redshift. We adopt a local, astrophysical
binary neutron star merger rate of 440 Gpc-3 y-1, the 90% confidence upper
bound from Ref. Abbott  _et al._ (2021)’s PowerLaw+Dip+Break model. The binary
tidal deformability for each event is determined from its masses by the
$m$–$\Lambda$ relation dictated by the equation of state. Neutron stars in the
twin-star mass range are randomly assigned to the hadronic or the hybrid
branch of the relation with equal probability.

We calculate optimal signal-to-noise ratios for the simulated mergers with
respect to three different detector networks: the “HLV” network includes LIGO-
Hanford, LIGO-Livingston and Virgo detectors at their design sensitivities;
the “A+” network upgrades the LIGO-Hanford and LIGO-Livingston detectors to A+
sensitivity; and the “XG” network replaces the LIGO detectors with Cosmic
Explorer detectors, and the Virgo detector with Einstein Telescope, at their
respective design sensitivities. For the HLV and A+ networks, we simulate
mergers within their projected luminosity distance ranges of 190 Mpc and 330
Mpc, respectively, for binary neutron stars Abbott  _et al._ (2018b). For the
XG network, whose binary neutron star range will extend to cosmological
distances Evans  _et al._ (2021); Borhanian and Sathyaprakash (2022), we
simulate only the nearby mergers ($z\lesssim 0.5$) that contribute the loudest
events. Our signal-to-noise ratio calculation is carried out with bilby Ashton
_et al._ (2019), using the IMRPhenomPv2_NRTidal waveform model Dietrich  _et
al._ (2019), a minimum frequency of 10 Hz (40 Hz for HLV), and a network
signal-to-noise ratio threshold of 12 for detection.

Given the assumed merger rate and distribution of sources, our simulated
populations yield expected astrophysical rates of 3, 15 and 5800 binary
neutron star mergers per year within the aforementioned ranges for HLV, A+ and
XG, respectively. The exact number of detections depends on the population
realization, but the detection efficiency is $\sim 60\%$ for A+ and $\sim
95\%$ for XG. Since the HLV network makes so few detections, we focus our
analysis on the A+ and XG scenarios. For each detected event, we simulate the
likelihoods in chirp mass $\mathcal{M}$, mass ratio $q$ and binary tidal
deformability $\tilde{\Lambda}$ as independent Gaussians, scaling their
standard deviations with the event’s signal-to-noise ratio in inverse
proportion to GW170817’s Farr  _et al._ (2016). We mock up the effect of the
detector noise realization by adding Gaussian noise that shifts the median of
the likelihood within one standard deviation.

An example of the resulting observed distribution of binary tidal
deformabilities vs chirp masses is illustrated in Fig. 2 for SKI5_2009. One
can clearly see the disjoint nature of the true, underlying distribution, but
it is less visible when the data is scattered by detector noise. For equations
of state with smaller $\Delta\Lambda$, the gaps in the distribution may not be
discernible by eye at all. The need to simultaneously account for statistical
uncertainty in the masses and tidal deformabilities of individual events, as
well as uncertainty in the underlying $m$–$\Lambda$ relation, motivates the
elaboration of the hierarchical inference framework we describe below.

##  IV Hierarchical Inference

We systematically search for evidence of twin stars in the binary tidal
deformability vs chirp mass distribution using a simple model. We prescribe a
parametric population model

$\pi(\mathcal{M},q,\tilde{\Lambda}|\boldsymbol{\theta},\boldsymbol{\lambda})=\pi(\tilde{\Lambda}|\mathcal{M},q;\boldsymbol{\theta})\pi(\mathcal{M},q|\boldsymbol{\lambda})$
(1)

built from a binary tidal deformability distribution

$\displaystyle\pi(\tilde{\Lambda}|\mathcal{M},q;\boldsymbol{\theta})=\int$
$\displaystyle
P(\tilde{\Lambda}|\Lambda_{1},\Lambda_{2},q)\pi(\Lambda_{1}|m_{1};\boldsymbol{\theta})\pi(\Lambda_{2}|m_{2};\boldsymbol{\theta})$
$\displaystyle\times
P(m_{1},m_{2}|\mathcal{M},q)d\Lambda_{1}d\Lambda_{2}dm_{1}dm_{2}$ (2)

that depends on parameters $\boldsymbol{\theta}$ describing an $m$–$\Lambda$
relation
$\pi(\Lambda|m;\boldsymbol{\theta})=\delta(\Lambda-\Lambda_{\boldsymbol{\theta}}(m))$,
and a mass distribution

$\pi(\mathcal{M},q;\boldsymbol{\lambda})=\int
P(\mathcal{M},q|m_{1},m_{2})\pi(m_{1},m_{2}|\boldsymbol{\lambda})dm_{1}dm_{2}$
(3)

that depends on parameters $\boldsymbol{\lambda}$, expressed here in terms of
a joint distribution $\pi(m_{1},m_{2}|\boldsymbol{\lambda})$ for the component
masses $m_{1,2}$. A hierarchical inference of the population model parameters
$\boldsymbol{\theta},\boldsymbol{\lambda}$ proceeds according to Mandel  _et
al._ (2019)

$\displaystyle
P(\boldsymbol{d}|\boldsymbol{\theta},\boldsymbol{\lambda})=\prod_{i}\frac{1}{\zeta(\boldsymbol{\theta},\boldsymbol{\lambda})}\int$
$\displaystyle P(d_{i}|\mathcal{M},q,\tilde{\Lambda})$
$\displaystyle\times\pi(\mathcal{M},q,\tilde{\Lambda}|\boldsymbol{\theta},\boldsymbol{\lambda})\,d\mathcal{M}\,dq\,d\tilde{\Lambda},$
(4)

where $P(d_{i}|\mathcal{M},q,\tilde{\Lambda})$ is the marginal gravitational-
wave likelihood in chirp mass, mass ratio and binary tidal deformability for
the $i$th observation, and $\zeta(\boldsymbol{\theta},\boldsymbol{\lambda})$
is the fraction of the population that is detected. Since we are only
interested in the $m$–$\Lambda$ relation parameters $\boldsymbol{\theta}$ in
our application, we fix the mass distribution as specified above. We also
ignore selection effects, as they do not impact the recovery of the parameters
of interest: in our population model, the detection fraction is approximately
independent of the equation of state.333The gravitational-wave selection
function for binary neutron star mergers is proportional to
$\mathcal{M}^{5/2}$ Chatziioannou and Farr (2020), so the detection fraction
depends on the equation of state through the maximum mass in the population,
set here by $M_{\rm TOV}$. However, $M_{\rm TOV}$ varies by only $\sim 10\%$
across the chosen equations of state. Thus, we approximate Eq. (IV) as

$P(\boldsymbol{d}|\boldsymbol{\theta})\propto\prod_{i}\int
P(d_{i}|\mathcal{M},q,\tilde{\Lambda})\,\pi(\mathcal{M},q,\tilde{\Lambda}|\boldsymbol{\theta})\,d\mathcal{M}\,dq\,d\tilde{\Lambda}.$
(5)

We adopt uniform priors on the $m$–$\Lambda$ relation parameters
$\boldsymbol{\theta}$, such that the posterior
$P(\boldsymbol{\theta}|\boldsymbol{d})$ is proportional to this equation.

To evaluate the likelihood $P(\boldsymbol{d}|\boldsymbol{\theta})$, we make a
Monte Carlo approximation for the integrals over $\mathcal{M}$ and $q$. The
remaining integral is resolved by the delta function in $\tilde{\Lambda}$ that
results from Eq. (IV). We then have

$\displaystyle P$ $\displaystyle(\boldsymbol{d}|\boldsymbol{\theta})\approx$
$\displaystyle\prod_{i}\sum_{j}P(d_{i}|\mathcal{M}_{j},q_{j},\tilde{\Lambda}_{\boldsymbol{\theta}}(\mathcal{M}_{j},q_{j}))\;|\;\mathcal{M}_{j},q_{j}\sim\pi(\mathcal{M},q)$
(6)

up to an overall normalization. We sample from this likelihood using a Markov-
chain Monte Carlo algorithm implemented with emcee Foreman-Mackey  _et al._
(2013).

To model the twin-star $m$–$\Lambda$ relations, we adopt the parametric form

$\Lambda_{\boldsymbol{\theta}}=\begin{cases}\Lambda_{1.4}^{\rm
had}\left(m/1.4\,M_{\odot}\right)^{-6}&1\,M_{\odot}\leq m\leq M_{t}\\\
\left[\Lambda_{1.4}^{\rm
had}\left(M_{t}/1.4\,M_{\odot}\right)^{-6}-\Delta\Lambda\right]\left[\left(\frac{m}{M_{t}}\right)^{-k}-\frac{(1-m/M_{t})}{(1-M_{\rm
TOV}/M_{t})}\left(\frac{M_{\rm TOV}}{M_{t}}\right)^{-k}\right]+\Lambda_{\rm
TOV}\frac{(1-m/M_{t})}{(1-M_{\rm TOV}/M_{t})}&M_{t}<m\leq M_{\rm TOV}\\\
\end{cases}$ (7)

with parameters $\boldsymbol{\theta}=\\{\Lambda_{1.4}^{\rm
had},M_{t},\Delta\Lambda,k,M_{\rm TOV},\Lambda_{\rm TOV}\\}$—the tidal
deformability of a canonical hadronic neutron star, the (median) twin-star
mass, the tidal deformability difference between twins of mass $M_{t}$, the
power-law slope of $\Lambda(m)$ on the hybrid branch, the maximum neutron star
mass and the tidal deformability of the maximum-mass star, respectively. To
keep the model simple, we do not model the finite extent of the twin star mass
range; $M_{t}$ simply represents its midpoint. Thus, the twin-star signature
that informs the hierarchical inference is the existence of a gap of width
$\Delta\Lambda$ in the tidal deformability distribution at $M_{t}$, rather
than the multimodality of the tidal deformability distribution itself. This
model for the mass-tidal deformability relation is shown in Fig. 3.

The uniform priors for the $m$–$\Lambda$ relation parameters are subjected to
the constraints $1.00\,M_{\odot}\leq M_{t}\leq 1.94\,M_{\odot}$,
$0\leq\Delta\Lambda\leq 1500$, $5\leq k\leq 12$. Finding that $\Lambda_{1.4}$
is always well-recovered, and that the results are relatively insensitive to
the maximum-mass parameters, we fix $M_{\rm TOV}$, $\Lambda_{1.4}$ and
$\Lambda_{\rm TOV}=3$ for simplicity.

Figure 3: Parametric model for the mass-tidal deformability relation used in
the hierarchical inference. The fit to the relation for the SKI5_2009 equation
of state is shown. The parameters we infer are highlighted in blue. The
parameters we fix for simplicity are shown in black.

##  V Results

We apply the hierarchical inference of twin-star parameters to the populations
simulated for the hybrid equations of state listed in Table 1, as well as the
purely hadronic equations of state SKI5 and SK272. For the hadronic cases, we
expect to recover the prior on $M_{t}$ and a posterior peaked at
$\Delta\Lambda=0$. For the hybrid cases supporting twin stars, we expect to
recover a posterior on $\Delta\Lambda$ peaked away from zero, and a constraint
on $M_{t}$. We consider the twin stars as detected if the marginal
$\Delta\Lambda$ highest-posterior-density interval excludes zero at 90%
confidence. As an alternative quantification of the evidence for twin stars,
we also compute an evidence ratio (Bayes factor) between the hypotheses that
twin stars are present and absent in the population.

The recovered twin star parameters for each equation of state are shown in
Fig. 4 for a sequence of observing scenarios, from one year at A+ sensitivity
to one month at XG sensitivity. We show the average recovery over 10
realizations of the population to mitigate statistical fluctuations in the
confidence regions. In the case of the equation of state with the strongest
phase transition, SKI5_2009, the one-dimensional marginal posterior on
$\Delta\Lambda$ already comfortably favors the presence of twin stars at 90%
confidence after one week of observation with the XG network. The
$\Delta\Lambda$ posterior for SKI5_2006, which has the same onset density but
a somewhat weaker phase transition, also satisfies our criterion for twin-star
detection after one week of XG observations.

For the equations of state with higher onset densities, and correspondingly
smaller $\Delta M_{t}$ and $\Delta\Lambda$, the tidal deformability resolution
in the XG observations is not sufficient to distinguish between zero and
finite $\Delta\Lambda$. However, after one month of XG observations,
$\Delta\Lambda$ larger than $\sim 100$ can be ruled out at 90% confidence for
both SK272_2506 and SK272_3505. This resolution limit in tidal deformability
also holds for the recoveries of the hadronic equations of state SKI5 and
SK272: the $\Delta\Lambda$ posteriors converge towards zero as expected, but
cannot exclude $\Delta\Lambda\lesssim 100$.

Hence, based on the 90% confidence intervals for $\Delta\Lambda$, we correctly
conclude after one month of XG observations that twin stars are present in the
binary neutron star populations for SKI5_2009 and SKI5_2006, and that twin
stars with $\Delta\Lambda\gtrsim 100$ are excluded in the other four
scenarios. In all cases, the $\Delta\Lambda$ posteriors after even two years
of A+ observations are essentially uninformative. This demonstrates that an XG
detector network is needed to definitively identify (or rule out) twin stars
with gravitational waves at the population level.

Figure 4: Posteriors on twin-star parameters of the mass-tidal deformability
relation for different equations of state. The 90% credible region of the
posterior, averaged over 10 population realizations, is shown for various
observing scenarios. The parameters corresponding to the injected equation of
state are indicated in black; in the bottom two panels, where the injected
equation of state is purely hadronic, the $M_{t}$ parameter has no intrinsic
meaning.

The evidence ratios between twin and no-twin hypotheses tell a similar story.
For each equation of state and observing scenario, we compute the Bayes factor
$B$ as a Savage-Dickey density ratio Dickey (1971), the ratio of marginal
posterior to marginal prior at $\Delta\Lambda=0$, obtaining the former via a
kernel density estimate. We plot the evolution of the Bayes factors in Fig. 5
as a function of the number of binary neutron star mergers detected within a
redshift $z<0.5$; the error bands quantify the variation in $B$ over 10
population realizations. As can be seen, for SKI5_2009 the evidence is
decidedly in favor of twin stars ($\log{B}<-6$) after a week of XG
observations ($\sim 100$ detections). For SKI5_2006, the evidence favors twin
stars more modestly ($\log{B}\approx-3$), although some population
realizations yield as decisive a result as for SKI5_2009. For the other four
cases, the Bayes factor test favors the no-twins hypothesis only moderately
($\log{B}\approx 3$) even after a month of XG observations ($\sim 500$
detections). This reinforces the conclusion that an XG network can confidently
detect $\Delta\Lambda\gtrsim 100$, but is resolution-limited for smaller tidal
deformability differences.

Figure 5: Model comparison between twin-star and no-twin hypotheses for each
injected equation of state. We show how the log Bayes factor evolves as a
function of the number of observations in the A+ era (light shading) and the
XG era (dark shading). A negative $\log B$ favors the presence of twin stars.
The shaded regions encompass 68% confidence intervals on the Bayes factor due
to its variation across 10 different population realizations.

Having established the prospects for twin-star detection, we now investigate
how well the twin-star parameters themselves can be measured. The main panels
of Fig. 4 relate the accuracy and precision with which $\Delta\Lambda$ and
$M_{t}$ are recovered. We find that the recoveries are accurate at the 90%
confidence level for all of the equations of state we study. However, for
SKI5_2009, the mode of the posterior slightly overestimates the tidal
deformability difference. We attribute this to the very small range $\Delta
M_{t}\approx 0.01\,M_{\odot}$ over which this equation of state supports twin
stars, which means that the gaps in the binary tidal deformability vs chirp
mass distribution are easily mimicked by statistical fluctuations in the data.
For SKI272_3505, the posterior mode instead slightly underestimates $M_{t}$,
which we attribute to the proximity between $M_{t}$ and $M_{\rm TOV}$ for this
equation of state. For the hadronic equations of state SKI5 and SK272, the
$M_{t}$ parameter is meaningless, and we observe that the recovered value
tends to be close to the midpoint of the neutron star mass spectrum.

The size of the 90% confidence regions in Fig. 4 demonstrates the increasing
precision in twin-star parameter measurements that can be expected from future
gravitational-wave observations. Overall, we find that $M_{t}$ is more
precisely constrained than $\Delta\Lambda$. In the best-case scenario, for
SKI5_2009, $M_{t}$ and $\Delta\Lambda$ are determined to within 1% and 15%
uncertainty, respectively, after one month of observation. In contrast, after
two years of A+ observations, the uncertainties are respectively 15% and 68%.
Even in the most pessimistic scenario, for SK272_3505, $M_{t}$ and
$\Delta\Lambda$ are measured with respective errors of 20% and 150% after one
month at XG sensitivity. Again, this illustrates the constraining power of XG
observations.

##  VI Discussion

Our investigations demonstrate how hierarchical inference on a population of
binary neutron star mergers can be used to systematically test for the
presence of twin stars. We find that an XG network may be able to discern twin
stars in the population after as little as one week of observation. Such a
network can resolve a tidal deformability difference between twins of several
hundred, and constrain the twin-star mass scale to within a few percent. When
realized in future, measurements of this kind will guide the development of
nuclear theory models for the neutron star equation of state.

While our analysis is not an exhaustive exploration of the parameter space for
twin stars, we forecast parameter constraints for a variety of equations of
state, including examples with first-order phase transitions of different
onset densities and strengths. Besides the choice of equation of state, our
forecasts also depend in principle on assumptions about the neutron star mass
distribution and the relative abundance of hadronic vs hybrid twins in the
detected population. We perform supplemental analyses in Appendix A that
indicate that these additional factors do not strongly impact the recovery of
the twin-star parameters.

Although the parametric $m$–$\Lambda$ relation model we adopt for the
hierarchical inference is deliberately kept simple, our results indicate that
it will suffice for accurate twin-star parameter measurements into the XG era.
Nonetheless, when future observing campaigns return thousands of binary
neutron star detections, it may be desirable to use a more sophisticated model
that does not gloss over the finite extent of the twin-star mass range. One
could also build the $m$–$\Lambda$ parameterization in terms of the
fundamental properties of the phase transition itself, such as the onset
density and the strength of the transition, in order to make measurements that
are more directly informative for nuclear theory. The general hierarchical
inference method laid out here can be adapted straightforwardly to accommodate
these modifications. It would also work to detect other phenomena giving rise
to ostensible twins, e.g. an overlap in the neutron star and black hole mass
spectra or a separate family of stellar-mass exotic compact objects.

###### Acknowledgements.

The authors thank Katerina Chatziioannou and Jolien Creighton for useful
suggestions about this work. P.L. is supported by the Natural Sciences &
Engineering Research Council of Canada (NSERC). K.C. is supported by the Czech
Academy of Sciences under the project number LQ100102101. The authors are
grateful for computational resources provided by the LIGO Lab and supported by
NSF Grants PHY-0757058 and PHY-0823459.

## References

  * Gerlach (1968) U. H. Gerlach, PhRv 172, 1325 (1968).
  * Kampfer (1981) B. Kampfer, JPhA 14, L471 (1981).
  * Glendenning and Kettner (2000) N. K. Glendenning and C. Kettner, A&A 353, L9 (2000), arXiv:astro-ph/9807155 [astro-ph] .
  * Schertler  _et al._ (2000) K. Schertler, C. Greiner, J. Schaffner-Bielich, and M. H. Thoma, NuPhA 677, 463 (2000), arXiv:astro-ph/0001467 [astro-ph] .
  * Alford  _et al._ (2013) M. G. Alford, S. Han, and M. Prakash, PhRvD 88, 083013 (2013), arXiv:1302.4732 [astro-ph.SR] .
  * Zdunik and Haensel (2013) J. L. Zdunik and P. Haensel, A&A 551, A61 (2013), arXiv:1211.1231 [astro-ph.SR] .
  * Benić  _et al._ (2015) S. Benić, D. Blaschke, D. E. Alvarez-Castillo, T. Fischer, and S. Typel, A&A 577, A40 (2015), arXiv:1411.2856 [astro-ph.HE] .
  * Alford  _et al._ (2015) M. G. Alford, G. F. Burgio, S. Han, G. Taranto, and D. Zappalà, PhRvD 92, 083002 (2015), arXiv:1501.07902 [nucl-th] .
  * Han and Steiner (2019) S. Han and A. W. Steiner, PhRvD 99, 083014 (2019), arXiv:1810.10967 [nucl-th] .
  * Schaffner-Bielich  _et al._ (2002) J. Schaffner-Bielich, M. Hanauske, H. Stöcker, and W. Greiner, PhRvL 89, 171101 (2002), arXiv:astro-ph/0005490 [astro-ph] .
  * Zacchi  _et al._ (2017) A. Zacchi, L. Tolos, and J. Schaffner-Bielich, PhRvD 95, 103008 (2017), arXiv:1612.06167 [astro-ph.HE] .
  * Alford and Sedrakian (2017) M. Alford and A. Sedrakian, PhRvL 119, 161104 (2017), arXiv:1706.01592 [astro-ph.HE] .
  * Christian  _et al._ (2018) J.-E. Christian, A. Zacchi, and J. Schaffner-Bielich, EPJA 54, 28 (2018), arXiv:1707.07524 [astro-ph.HE] .
  * Oppenheimer and Volkoff (1939) J. R. Oppenheimer and G. M. Volkoff, PhRv 55, 374 (1939).
  * Tolman (1939) R. C. Tolman, PhRv 55, 364 (1939).
  * Flanagan and Hinderer (2008) É. É. Flanagan and T. Hinderer, PhRvD 77, 021502 (2008), arXiv:0709.1915 [astro-ph] .
  * Hinderer (2008) T. Hinderer, ApJ 677, 1216 (2008), arXiv:0711.2420 [astro-ph] .
  * Landry and Poisson (2014) P. Landry and E. Poisson, PhRvD 89, 124011 (2014), arXiv:1404.6798 [gr-qc] .
  * Paschalidis  _et al._ (2018) V. Paschalidis, K. Yagi, D. Alvarez-Castillo, D. B. Blaschke, and A. Sedrakian, PhRvD 97, 084038 (2018), arXiv:1712.00451 [astro-ph.HE] .
  * Montaña  _et al._ (2019) G. Montaña, L. Tolós, M. Hanauske, and L. Rezzolla, PhRvD 99, 103009 (2019), arXiv:1811.10929 [astro-ph.HE] .
  * Christian  _et al._ (2019) J.-E. Christian, A. Zacchi, and J. Schaffner-Bielich, PhRvD 99, 023009 (2019), arXiv:1809.03333 [astro-ph.HE] .
  * Pang  _et al._ (2020) P. T. H. Pang, T. Dietrich, I. Tews, and C. Van Den Broeck, PhRvR 2, 033514 (2020), arXiv:2006.14936 [astro-ph.HE] .
  * Wang  _et al._ (2022) Q.-W. Wang, C. Shi, Y. Yan, and H.-S. Zong, NuPhA 1025, 122489 (2022), arXiv:1912.02312 [astro-ph.HE] .
  * Abbott  _et al._ (2017) B. P. Abbott, R. Abbott, T. D. Abbott, F. Acernese, K. Ackley, _et al._ , Phys. Rev. Lett. 119, 161101 (2017), arXiv:1710.05832 [gr-qc] .
  * Abbott  _et al._ (2019) B. P. Abbott, R. Abbott, T. D. Abbott, F. Acernese, K. Ackley, _et al._ , PhRvX 9, 011001 (2019), arXiv:1805.11579 [gr-qc] .
  * Abbott  _et al._ (2018a) B. P. Abbott, R. Abbott, T. D. Abbott, F. Acernese, K. Ackley, _et al._ , Phys. Rev. Lett. 121, 161101 (2018a), arXiv:1805.11581 [gr-qc] .
  * Antoniadis  _et al._ (2013) J. Antoniadis, P. C. C. Freire, N. Wex, T. M. Tauris, R. S. Lynch, _et al._ , Sci 340, 448 (2013), arXiv:1304.6875 [astro-ph.HE] .
  * Cromartie  _et al._ (2020) H. T. Cromartie, E. Fonseca, S. M. Ransom, P. B. Demorest, Z. Arzoumanian, _et al._ , NatAs 4, 72 (2020), arXiv:1904.06759 [astro-ph.HE] .
  * Fonseca  _et al._ (2021) E. Fonseca, H. T. Cromartie, T. T. Pennucci, P. S. Ray, A. Y. Kirichenko, _et al._ , ApJL 915, L12 (2021), arXiv:2104.00880 [astro-ph.HE] .
  * Christian and Schaffner-Bielich (2020) J.-E. Christian and J. Schaffner-Bielich, ApJL 894, L8 (2020), arXiv:1912.09809 [astro-ph.HE] .
  * Riley  _et al._ (2019) T. E. Riley, A. L. Watts, S. Bogdanov, P. S. Ray, R. M. Ludlam, _et al._ , ApJL 887, L21 (2019), arXiv:1912.05702 [astro-ph.HE] .
  * Miller  _et al._ (2019) M. C. Miller, F. K. Lamb, A. J. Dittmann, S. Bogdanov, Z. Arzoumanian, _et al._ , ApJL 887, L24 (2019), arXiv:1912.05705 [astro-ph.HE] .
  * Christian and Schaffner-Bielich (2022) J.-E. Christian and J. Schaffner-Bielich, ApJ 935, 122 (2022), arXiv:2109.04191 [astro-ph.HE] .
  * Riley  _et al._ (2021) T. E. Riley, A. L. Watts, P. S. Ray, S. Bogdanov, S. Guillot, _et al._ , ApJL 918, L27 (2021), arXiv:2105.06980 [astro-ph.HE] .
  * Miller  _et al._ (2021) M. C. Miller, F. K. Lamb, A. J. Dittmann, S. Bogdanov, Z. Arzoumanian, _et al._ , ApJL 918, L28 (2021), arXiv:2105.06979 [astro-ph.HE] .
  * Li  _et al._ (2021) J. J. Li, A. Sedrakian, and M. Alford, PhRvD 104, L121302 (2021), arXiv:2108.13071 [astro-ph.HE] .
  * Adhikari  _et al._ (2021) D. Adhikari, H. Albataineh, D. Androic, K. Aniol, D. S. Armstrong, _et al._ , PhRvL 126, 172502 (2021), arXiv:2102.10767 [nucl-ex] .
  * Christian and Schaffner-Bielich (2021) J.-E. Christian and J. Schaffner-Bielich, PhRvD 103, 063042 (2021), arXiv:2011.01001 [astro-ph.HE] .
  * Tsaloukidis  _et al._ (2022) L. Tsaloukidis, P. S. Koliogiannis, A. Kanakis-Pegios, and C. C. Moustakidis, arXiv , arXiv:2210.15644 (2022), arXiv:2210.15644 [astro-ph.HE] .
  * Aasi  _et al._ (2015) J. Aasi, B. P. Abbott, R. Abbott, T. Abbott, M. R. Abernathy, _et al._ , Classical and Quantum Gravity 32, 074001 (2015), arXiv:1411.4547 [gr-qc] .
  * Acernese  _et al._ (2015) F. Acernese, M. Agathos, K. Agatsuma, D. Aisa, N. Allemandou, _et al._ , Classical and Quantum Gravity 32, 024001 (2015), arXiv:1408.3978 [gr-qc] .
  * Akutsu  _et al._ (2021) T. Akutsu, M. Ando, K. Arai, Y. Arai, S. Araki, _et al._ , PTEP 2021, 05A102 (2021), arXiv:2009.09305 [gr-qc] .
  * Abbott  _et al._ (2018b) B. P. Abbott, R. Abbott, T. D. Abbott, M. R. Abernathy, F. Acernese, _et al._ , LRR 21, 3 (2018b), arXiv:1304.0670 [gr-qc] .
  * Colombo  _et al._ (2022) A. Colombo, O. S. Salafia, F. Gabrielli, G. Ghirlanda, B. Giacomazzo, _et al._ , ApJ 937, 79 (2022), arXiv:2204.07592 [astro-ph.HE] .
  * Patricelli  _et al._ (2022) B. Patricelli, M. G. Bernardini, M. Mapelli, P. D’Avanzo, F. Santoliquido, _et al._ , MNRAS 513, 4159 (2022), arXiv:2204.12504 [astro-ph.HE] .
  * Evans  _et al._ (2021) M. Evans, R. X. Adhikari, C. Afle, S. W. Ballmer, S. Biscoveanu, _et al._ , arXiv , arXiv:2109.09882 (2021), arXiv:2109.09882 [astro-ph.IM] .
  * Maggiore  _et al._ (2020) M. Maggiore, C. Van Den Broeck, N. Bartolo, E. Belgacem, D. Bertacca, _et al._ , JCAP 2020, 050 (2020), arXiv:1912.02622 [astro-ph.CO] .
  * Borhanian and Sathyaprakash (2022) S. Borhanian and B. S. Sathyaprakash, arXiv , arXiv:2202.11048 (2022), arXiv:2202.11048 [gr-qc] .
  * Smith  _et al._ (2021) R. Smith, S. Borhanian, B. Sathyaprakash, F. Hernandez Vivanco, S. E. Field, _et al._ , PhRvL 127, 081102 (2021), arXiv:2103.12274 [gr-qc] .
  * Chatziioannou and Han (2020) K. Chatziioannou and S. Han, PhRvD 101, 044019 (2020), arXiv:1911.07091 [gr-qc] .
  * Agrawal  _et al._ (2003) B. K. Agrawal, S. Shlomo, and V. Kim Au, PhRvC 68, 031304 (2003), arXiv:nucl-th/0308042 [nucl-th] .
  * Reinhard and Flocard (1995) P. G. Reinhard and H. Flocard, NuPhA 584, 467 (1995).
  * Bertulani and Valencia (2019) C. A. Bertulani and J. Valencia, PhRvC 100, 015802 (2019), arXiv:1904.01078 [nucl-th] .
  * Reed  _et al._ (2021) B. T. Reed, F. J. Fattoyev, C. J. Horowitz, and J. Piekarewicz, PhRvL 126, 172503 (2021), arXiv:2101.03193 [nucl-th] .
  * Landry and Read (2021) P. Landry and J. S. Read, ApJL 921, L25 (2021), arXiv:2107.04559 [astro-ph.HE] .
  * Abbott  _et al._ (2021) R. Abbott, T. D. Abbott, F. Acernese, K. Ackley, C. Adams, _et al._ , arXiv , arXiv:2111.03634 (2021), arXiv:2111.03634 [astro-ph.HE] .
  * Madau and Dickinson (2014) P. Madau and M. Dickinson, ARA&A 52, 415 (2014), arXiv:1403.0007 [astro-ph.CO] .
  * Ashton  _et al._ (2019) G. Ashton, M. Hübner, P. D. Lasky, C. Talbot, K. Ackley, _et al._ , ApJS 241, 27 (2019), arXiv:1811.02042 [astro-ph.IM] .
  * Dietrich  _et al._ (2019) T. Dietrich, S. Khan, R. Dudi, S. J. Kapadia, P. Kumar, _et al._ , PhRvD 99, 024029 (2019), arXiv:1804.02235 [gr-qc] .
  * Farr  _et al._ (2016) B. Farr, C. P. L. Berry, W. M. Farr, C.-J. Haster, H. Middleton, _et al._ , ApJ 825, 116 (2016), arXiv:1508.05336 [astro-ph.HE] .
  * Mandel  _et al._ (2019) I. Mandel, W. M. Farr, and J. R. Gair, Mon. Not. Royal Astron. Soc. 486, 1086 (2019), arXiv:1809.02063 [physics.data-an] .
  * Chatziioannou and Farr (2020) K. Chatziioannou and W. M. Farr, PhRvD 102, 064063 (2020), arXiv:2005.00482 [astro-ph.HE] .
  * Foreman-Mackey  _et al._ (2013) D. Foreman-Mackey, D. W. Hogg, D. Lang, and J. Goodman, PASP 125, 306 (2013), arXiv:1202.3665 [astro-ph.IM] .
  * Dickey (1971) J. M. Dickey, The Annals of Mathematical Statistics 42, 204 (1971).
  * Farr and Chatziioannou (2020) W. M. Farr and K. Chatziioannou, Research Notes of the American Astronomical Society 4, 65 (2020), arXiv:2005.00032 [astro-ph.GA] .

##  Appendix A Supplemental Material

In this appendix, we revisit some of the choices made in the main analysis and
demonstrate that our conclusions are robust against alternative assumptions.
We perform posterior predictive checks to reinforce the accuracy of our
inferences. We also examine the connection between the $m$–$\Lambda$ relation
for hybrid stars and the morphology of the binary tidal deformability vs chirp
mass distribution in greater detail.

###  A.1 Morphology of the tidal deformability distribution

Figure 6: Twin-star–supporting equations of state used in this work. The
discontinuities in baryon density at fixed pressure are first-order phase
transitions of Maxwell type. The phase transitions occurring in these
equations of state are strong enough to give rise to stable twin stars.

In the main analysis, we argued that an equation of state’s support for twin
stars translates into a disjoint two-dimensional distribution of binary tidal
deformability

$\tilde{\Lambda}=\frac{16}{13}\frac{(1+12q)\Lambda_{1}+(q+12)q^{4}\Lambda_{2}}{(1+q)^{5}}$
(8)

vs chirp mass

$\mathcal{M}=\frac{(m_{1}m_{2})^{3/5}}{(m_{1}+m_{2})^{1/5}},$ (9)

whereas the equivalent distribution is contiguous for an equation of state
without a disconnected hybrid star branch. The gaps in the distribution
originate directly from the unstable segment of the $m$–$\Lambda$ relation
that connects the hadronic and hybrid branches for twin-star-supporting
equations of state: these are the dashed segments in Fig. 1. They correspond
to strong first-order phase transitions in the equation of state, where the
baryon density jumps discontinuously at fixed pressure—see Fig. 6. The gaps in
the $\tilde{\Lambda}$ vs $\mathcal{M}$ distribution occur where, along lines
of constant mass ratio $q=m_{2}/m_{1}$, one component mass lies on the
unstable segment. This is illustrated in Fig. 7, which breaks down the
distribution in Fig. 2 along lines of fixed mass ratio. The upper
(respectively, lower) gap in the distribution for SKI5_2009 corresponds to
$m_{2}$ ($m_{1}$) lying on the unstable segment. Note how a similar breakdown
of the distribution for the purely hadronic SKI5 equation of state does not
evince the appearance of any gaps.

Figure 7 also shows the same distribution for SK272_3505 and its purely
hadronic counterpart, SK272. In this case, due to its very small $\Delta
M_{t}$ and $\Delta\Lambda$, the gaps in the distribution for SK272_3505 are
barely visible at scale. However, along lines of constant mass ratio, one can
see the distribution is broken into distinct segments. Again, this contrasts
with the continuous distribution for SK272.

---
Figure 7: The joint binary tidal deformability and chirp mass distribution
probed by one month of XG observations for selected equations of state studied
in this work. The distribution is broken down according to mass ratio, showing
how the gaps that indicate the presence of twin stars arise. A uniform neutron
star mass distribution is assumed.

###  A.2 Posterior predictive checks

|
---|---
|
Figure 8: Posterior predictive checks for selected twin-star-supporting
equations of state. We compare the posterior predictive distribution (PPD) for
the binary tidal deformability vs chirp mass distribution (left) and the
mass–tidal deformability relation (right) to the injected values. The
posterior predictive distribution is constructed from the approximate mode of
the recovered three-dimensional posterior in $\Delta\Lambda$, $M_{t}$ and $k$.

To bolster the results of the main analysis, we perform posterior predictive
checks and examine how constraints on the twin-star parameters evolve with the
number of binary neutron star detections. For each injected equation of state,
we extract the approximate maximum-likelihood parameters
($\Delta\Lambda$,$M_{t}$,$k$) from the population-averaged posterior after one
month of XG observations. We use these parameters to reconstruct the best-fit
$m$–$\Lambda$ relation according to Eq. (7). We then sample a uniform binary
neutron star population and reconstruct the predicted distribution of binary
tidal deformability vs chirp mass for the population. Figure 8 compares these
predictions for the $m$–$\Lambda$ relation and the $\tilde{\Lambda}$ vs
$\mathcal{M}$ distribution to their true values for two selected equations of
state, SKI5_2009 and SK272_3505.

For SKI5_2009, we observe that the predicted distributions are close, but not
perfect, matches to the injected ones. The main discrepancy occurs at large
masses, where the model’s fixed power-law slope does not allow for the
flexibility required to track the actual $m$ vs $\Lambda$ trend. This is
essentially a deliberate tradeoff in the model, which prizes simplicity over
fidelity. Similarly, by construction, the model does not track the
$m$–$\Lambda$ relation right through the unstable segment connecting the
hadronic and hybrid branches. However, our simplified treatment of this
juncture manifestly reproduces the right $\tilde{\Lambda}$ vs $\mathcal{M}$
morphology and accurately locates the discontinuity in the $m$–$\Lambda$
relation.

For SK272_3505, which has a much higher phase transition onset density and
much smaller $\Delta M_{t}$ and $\Delta\Lambda$, the model clearly struggles
to match both the location and extent of the hybrid branch. The mode of the
posterior, upon which the illustrated posterior predictive check is based,
favors a good fit to the latter at the expense of the former. However, as this
equation of state is the one with the largest uncertainties a posteriori, the
posterior encompasses other parameter combinations that recover the twin star
mass scale more accurately.

This can be seen in the evolution of the twin star parameter uncertainties as
a function of the number of $z<0.5$ binary neutron star detections. In Fig. 9,
we show the posterior 90% confidence intervals on $\Delta\Lambda$ and $M_{t}$
for the scenarios that are presented in Fig. 4. The parameters are recovered
accurately for all of the equations of state, except for a slight
underestimate of $M_{t}$ in the case of SK272_3505, which is related to the
mismatch in the posterior predictive check. The precision in the parameter
measurements increases significantly when passing from the A+ to the XG
detector network.

|
---|---
Figure 9: Constraints on the tidal deformability difference between twins and
the twin-star mass scale as a function of the number of binary neutron star
mergers detected within a redshift $z<0.5$. We show the evolution of the 90%
confidence contours. The era of A+ (respectively, XG) observations is shown
with lighter (darker) shading. We average over 10 population realizations to
mitigate fluctuations in statistical uncertainties.

###  A.3 Effect of the population realization

Figure 10: Posteriors on twin-star parameters of the mass-tidal deformability
relation for the SKI5_2009 equation of state. The 90% credible region of the
posterior is shown for the best-case A+ and XG observing scenarios we
consider. We show how the credible regions differ across several different
population realizations (faint traces) compared to the population averages
(solid lines) that were displayed in the equivalent panel of Fig. 4. Note how
the variation due to the population realization diminishes as more
observations are accumulated.

The results presented in the main analysis average over 10 realizations of the
binary neutron star population. Here we show how much the parameter
constraints vary with the population realization. In Fig. 10, we show the two-
dimensional posterior 90% confidence contours on $\Delta\Lambda$ and $M_{t}$
for SKI5_2009 for five population realizations after two years of A+
observations and one month of XG observations—these are the best-case A+ and
XG observing scenarios, respectively. One can see that statistical
fluctuations due to the population realization can significantly impact the A+
contours; however, as the detector network improves and the number of
observations increases, this variability becomes less significant.

###  A.4 Effect of the branching ratio

Figure 11: Posteriors on twin-star parameters of the mass-tidal deformability
relation for the SKI5_2009 equation of state, under a different assumption
about the hadronic vs hybrid branching ratio for twin stars. The 90% credible
region of the posterior, averaged over 10 population realizations, is shown
for various observing scenarios. Note the qualitative similarity with the
equivalent panel of Fig. 4.

We assumed above that neutron stars with masses in the twin star range lie
with equal probability on the hadronic or hybrid branch of the $m$–$\Lambda$
relation. Here we make an alternative choice and show that the twin-star
parameter constraints are largely unchanged in our best-case scenario. We
assume that twin stars are distributed uniformly in central density, such that
the probability for an arbitrary twin star to lie on the hadronic vs the
hybrid branch is proportional to the ratio of density ranges supported by each
twin-star branch. For instance, if hadronic twins span a central density range
of $0.3\,\rho_{\rm nuc}$, and hybrid twins span only $0.1\,\rho_{\rm nuc}$,
there should be three times as many hadronic twins as hybrid ones. Repeating
the main analysis for the SKI5_2009 equation of state under these conditions,
we obtain Fig. 12’s population-averaged posterior on the parameters
$\Delta\Lambda$ and $M_{t}$. As was the case under the original branching
ratio assumption, twin stars are definitively identified in the population
after one week of XG observations.

###  A.5 Effect of the mass distribution

Figure 12: Posteriors on twin-star parameters of the mass-tidal deformability
relation for different equations of state, under the assumption of a bimodal
neutron star mass distribution. The 90% credible region of the posterior,
averaged over 10 population realizations, is shown for various observing
scenarios. Note the qualitative similarity with Fig. 4.

In the main analysis, we adopted a uniform mass distribution for the binary
neutron star population. We now revisit that assumption and demonstrate that
our conclusions are unchanged if we instead use the bimodal neutron star mass
distribution from Ref. Farr and Chatziioannou (2020), inspired by observations
of Galactic pulsars. We continue to assume that neutron stars pair randomly
into binaries. When we repeat the main analysis with this new population
model, we obtain the twin-star parameter constraints shown in Fig. 12. As
above, the results are averaged over 10 population realizations. We obtain
qualitatively similar constraints as with the uniform neutron star mass
distribution: one month of XG observations can identify the presence of twin
stars in the SKI5_2009 and SKI5_2006 scenarios, while $\Delta\Lambda\gtrsim
100$ is ruled out in the other four scenarios.