# Knowledge-based in silico models and dataset for the comparative evaluation
of mammography AI for a range of breast characteristics, lesion conspicuities
and doses

E. Sizikova, N. Saharkhiz, D. Sharma, M. Lago, B. Sahiner, J. G. Delfino, A.
Badano
Office of Science and Engineering Laboratories
Center for Devices and Radiological Health
U.S. Food and Drug Administration
Silver Spring, MD 20993 USA

###### Abstract

To generate evidence regarding the safety and efficacy of artificial
intelligence (AI) enabled medical devices, AI models need to be evaluated on a
diverse population of patient cases, some of which may not be readily
available. We propose an evaluation approach for testing medical imaging AI
models that relies on in silico imaging pipelines in which stochastic digital
models of human anatomy (in object space) with and without pathology are
imaged using a digital replica imaging acquisition system to generate
realistic synthetic image datasets. Here, we release M-SYNTH***Code and data
links available at: https://github.com/DIDSR/msynth-release/, a dataset of
cohorts with four breast fibroglandular density distributions imaged at
different exposure levels using Monte Carlo x-ray simulations with the
publicly available Virtual Imaging Clinical Trial for Regulatory Evaluation
(VICTRE) toolkit. We utilize the synthetic dataset to analyze AI model
performance and find that model performance decreases with increasing breast
density and increases with higher mass density, as expected. As exposure
levels decrease, AI model performance drops with the highest performance
achieved at exposure levels lower than the nominal recommended dose for the
breast type.

Figure 1: Overview of the computational pipeline components for generating the
M-SYNTH in silico dataset for medical imaging AI evaluation.

##  1 Introduction

The goal of this work is to demonstrate that AI models for medical imaging can
be evaluated using simulations, specifically, using an in silico (also known
as synthetic) imaging pipeline equipped with a stochastic model for human
anatomy and disease [1]. We show that in silico methods can constitute rich
sources of data with realistic physical variability for performing comparative
analysis of AI device performance.

To date, computational models have been applied to some extent for the
analysis of nearly all medical imaging modalities and for a wide variety of
clinical tasks [2]. Since it is critical to ensure patient safety and system
effectiveness in healthcare applications, rigorous and thorough testing
procedures must be performed in order to study performance in the intended
population including subpopulations of interest. To prevent estimates that
might be biased by overfitting, model testing is typically performed on a
previously unseen dataset. However, datasets consisting of patient images may
present a limited distribution of the variability in human anatomy and may not
always capture rare, but life-critical cases, and may be biased towards
specific populations or parameters of image acquisition devices dominant at
specific clinical sites. In addition, patient data and associated health
records may not be available due to patient privacy, cost, or additional risk
associated with additional imaging procedures. Precise mass location and
extent (e.g., mass boundaries) are typically not available in the patient’s
records, and it is burdensome, error-prone, and sometimes impossible to
collect this information retrospectively. In many medical imaging
applications, these limitations pose a significant barrier to development and
evaluation of novel computational techniques in medical imaging products.

We propose evaluating AI models using physics-based simulations. We create
realistic test cases by imaging digital objects using digital image
acquisition systems. Our in silico testing pipeline offers the ability to
control both object and acquisition parameters, and generate highly realistic
test cases (see Figure 1). We show that digital objects and computer simulated
replicas of image acquisition devices offer a rich source of realistic data
capturing a variety of patient and imaging conditions for evaluation purposes.
In particular, our approach (and associated dataset) allows for performing
comparative analysis of AI performance across physical breast properties
(e.g., mass size) and imaging characteristics (e.g., radiation dose). Such
testing typically cannot be performed with patient data, as the data may be
too costly to collect or unsafe to acquire (e.g., one cannot ethically re-
image the same patient multiple times using ionizing radiation). Our
contributions in this work can be summarized as follows:

  * •

We demonstrate that, using this approach, we can detect differences in AI
model performance based on selected image acquisition device or physical
object model parameters. Specifically, we evaluate the effect of image
acquisition (radiation dose) and object model (breast and mass densities, mass
size) parameters on the performance of the AI model.

  * •

We release a dataset, M-SYNTH, to facilitate testing with pre-computed data
using the proposed pipeline. The dataset consists of 1,200 stochastic
knowledge-based models and their associated digital mammography (DM) images
with varying physical (breast density, mass size and density) and imaging
(dose) characteristics.

##  2 Background

First, we introduce the concepts of knowledge-based models and physics-based
imaging simulation that form the in silico imaging pipeline, the foundation of
our work.

Object Models. Knowledge-based (KB) models incorporate information about the
physical world into the data generation process to create realistic virtual
representations of human parts or organs [3]. As discussed in [1], large
cohorts of digital stochastic human models can be represented by:

$\\{\mathbf{f}_{s}\\}_{s=1}^{S}=\sum_{n}\theta_{n}^{s}\phi_{n}(\boldsymbol{r}),$
(1)

where $s$ denotes a particular state or random realization of a digital human
in a cohort of size $S$, $\boldsymbol{r}$ denotes a spatial variable,
$\phi_{n}$ denote expansion (basis) functions, and $\theta_{n}$ denote
expansion coefficients. Knowledge-based models specifically are constructed by
sampling a set of $\theta_{n}$ in Eq. 1 from distributions representing the
relevant model characteristics, given a specific $\phi_{n}$ based on the
application. The characteristics of the distributions are often derived from
physical or biological measurements. In the case of breast, knowledge-based
models allow us to vary physical patient characteristics including breast
size, breast shape, mass size and mass density (see Figures 4, 4 and 4).

Specifically, the object (breast) is a model $D$, parameterized by a vector
$x$ characterizing a fixed, user-defined set of physiological properties
(e.g., breast density, mass presence, mass size, glandularity). Given a sample
$x_{s}$, we can generate a realistic, high-resolution object $f_{s}=D(x_{s})$.
We rely on Graff’s breast model [3] as the KB model for this project and
describe its properties in Section 3.

Digital Mammography (DM) image generation. Once created, KB models are imaged
using simulations of x-ray transport through the materials present in each KB
model. The image acquisition device $I$ is a parametric model that receives
the object $d_{i}$ as well as user-defined choices for control parameters $y$
(e.g., detector type, radiation dose) and outputs an image
$r_{i,j}=I(d_{i},y_{j})$ given a sample choice of parameters $y_{j}$ and an
input object $d_{i}$. Parameters of such a system (e.g., geometry, source
characteristics, detector technology, anti-scatter grid, etc.) can emulate
system geometries and x-ray acquisition parameters found in commercially
available imaging device (e.g., mammography) specifications. In our work, we
used MC-GPU [4], a Monte Carlo x-ray simulation software implemented on GPUs
that generates mammography images. Additional details for this component of
the pipeline can be found in Section 4.

Related work in generative image models. The in silico imaging pipeline
described above is highly related to medical imaging generation using
generative models. One popular type of generative model is a generative
adversarial network (GAN) [5], which learns a mapping from a low-dimensional
representation to images at resolution. Generative models have been applied to
a variety of medical image generation tasks [6]. For example, Guan [7] showed
that GAN-generated synthetic images can be used to augment a smaller patient
breast image dataset for breast image classification. [8] introduced image-
based GAN to generate high resolution images conditioned on pixel-level mask
constraints. GANs may not correctly capture the link between input parameters
and outputs, and thus, are prone to generating unrealistic examples [9]. A
number of alternative types of generative models [10, 11, 12, 13, 14] have
been developed that address its limitations, such as training instabilities
and unrealistic output images. A key advantage of generative models is that
their run time can be faster than fully-detailed, object-space simulations,
and it remains important to explore and compare both techniques. Their key
limitation is that they require large training datasets and typically learn
noise and artifacts from the imaging system [15]. In particular, all image
acquisition systems have a null space, i.e., the set of object-space details
that are not observed in the acquired images due to imaging system limitations
(e.g., finite spatial and temporal resolution). Null space constraints limit
the ability of generative models to describe certain components of patient
anatomy and pathology. Simulation-based testing has been proposed in other
fields, such as autonomous vehicle navigation [16], and is related to the
concept of generating adversarial perturbations in the image [17, 18, 19] and
the physical property space [20, 21, 22]. For example, [23] introduced 3DB, a
photo-realistic simulation framework to debug and improve computer vision
models. Inspired by these works, we propose to evaluate medical imaging AI
using images generated using KB models and physics simulations and release a
dataset to facilitate such exploration.

##  3 Dataset Generation

The use of in silico imaging allows for the generation of large object and
image datasets without the need of human clinical trials. Here, we take
advantage of the benefits of the in silico approach to perform comparative
analysis of AI model performance across different physical properties of the
case population of breast models. We rely on the VICTRE pipeline †††See VICTRE
Github Page and FDA Regulatory Science Tools (RST) Catalog. for generating
breast models and their corresponding DM images. Previous work [24] has shown
that the VICTRE pipeline replicated the results of a clinical study comparing
DM and digital breast tomosynthesis (DBT) involving hundreds of enrolled
women. An overview of the data generation process can be seen in Figure 1.

Breast Model Synthesis. In silico breast models [3] (also known as breast
imaging phantoms) were generated using a procedural analytic model which
allows for adjusting various patient characteristics including breast shape,
size and glandular density. The models are compressed in the craniocaudal
direction using FeBio [25], an open source finite-element software. We
simplified the breast materials in non-glandular (as fat) or glandular tissue
with Young’s modulus and Poisson ratio of $E=5Pa$, $\nu=0.49$ and $E=15Pa$,
$\nu=0.49$, respectively. Lesions were inserted in a subset to create the
signal-present cohort. These models were then imaged using a state-of-the-art
Monte Carlo x-ray transport code (MC-GPU) [4].

We studied breast densities of extremely dense (referred to as “dense”),
heterogeneously dense (referred to as “hetero”), scattered, and fatty,
matching the distributions from [24]. For each breast density, a different
breast size is used to correspond with population statistics. Therefore, the
dense breast is the smallest, followed by heterogeneously dense, then
scattered, and then fatty. Each breast model was compressed to 3.5 cm, 4.5 cm,
5.5 cm, and 6.0 cm for each respective density, mimicking the organ
compression during the imaging. Random spiculated breast masses were generated
using the de Sisternes model [26] with three different sizes (5 mm, 7 mm and 9
mm radii) and mass density was set to be a factor of glandular tissue density
(1.0, 1.06 and 1.1 times). Note that for dense and hetero breasts, we only
used mass sizes of 5 and 7 mm, since 9 mm masses do not fit within the breast
region. No micro-calcification clusters were inserted. To create the signal-
present cohort, a single spiculated mass was inserted in half of the cases at
randomly chosen locations chosen from a list of candidate sites determined by
the position of the terminal duct lobular units. The resulting in silico
dataset comprises of 1,200 digital breast models, corresponding to 300
patients per breast size/density. Compared to the original VICTRE trial [24],
we introduce variations in mass size and density. Samples of model
realizations are shown in Figures 4, 4 and 4. Note that the bounding boxes are
only to make the masses more conspicuous for visualization purposes only.

(a) 5 mm (b) 7 mm (c) 9 mm (d) 5 mm (e) 7 mm (f) 9 mm

Figure 2: Effect of varying mass size (5 mm to 9 mm radius) in a fatty breast.
Two breast models are shown, first: (a)-(c), and second: (d)-(f). Dose (# of
hist.) $2.22\times 10^{10}$ and mass density 1.1 remain constant. Bounding
boxes are placed here to indicate the location of the masses.

(a) 1.0 (b) 1.06 (c) 1.1 (d) 1.0 (e) 1.06 (f) 1.1

Figure 3: Effect of varying mass density (1.0 to 1.1 times glandular tissue
density) in a fatty breast. Two models are shown, first: (a)-(c), and second:
(d)-(f). Dose (# of hist.) $2.22\times 10^{10}$ and mass size 7 mm remain
constant. Bounding boxes are placed here to indicate the location of the
masses.

(a) Fatty (b) Scattered (c) Hetero

(d) Dense (e) Renderings

Figure 4: Cohort variability. Varying breast density: (L to R) Fatty,
Scattered, Heterogeneously dense, and Dense with mass size 7 mm and mass
density 1.1. Note that dose changes with breast density. (e) shows artistic
renderings of models for each composition (details in Kim et al. [27].

Digital Mammography (DM) Generation. To simulate the x-ray imaging process, we
used MC-GPU [4], a Monte Carlo x-ray simulation software implemented on GPUs
that generates DM images. The detector model relies on system geometries and
x-ray acquisition parameters inspired by the currently available Siemens
Mammomat Inspiration DM system. The dosimetric and x-ray acquisition
parameters were selected based on publicly available device specifications and
clinical recommendations for each compressed breast thickness and
glandularity. We applied 20-100% of the clinically recommended dose for each
breast density. See Badal et al. [4] for the exact parameter values and doses
delivered to each breast and Sengupta et al. [28] for additional details.
X-ray photons arriving at the detector are tracked until first photoelectric
interaction incorporating fluorescence effects by generating and tracking a
secondary x-ray based on the fluorescence yield in a uniformly random
direction. Electronic noise is added to the pixel variance. The focal spot
blurring in the source was modeled as a 3D Gaussian probability distribution
with a full-width-at-half-maximum of 300 $\mu$m. A tungsten anode filtered
with 50 $\mu$m rhodium was used with a peak voltage of 28 kV for fatty and
scattered breasts and 30 kV for dense and heterogeneously dense breasts. The
same analytical anti-scatter grid was also included for generating the DM
images. (5:1 ratio, 31 line pairs/mm), see [4]. The resulting detector model
(known as DIR in [28]) is representative of a solid-state amorphous selenium
transducer in a direct detector configuration. Visualizations of generated
images and masses can be seen in Figure 5. A summary of complete parameters
used to generate data points in the presented dataset is described in Table 1.
In Figure 7, we report statistics of dose levels corresponding to the dataset.

(a) $4.44\times 10^{9}$ (b) $8.88\times 10^{9}$ (c) $1.33\times 10^{10}$ (d)
$1.78\times 10^{10}$ (e) $2.22\times 10^{10}$

Figure 5: Imaging. Effect of increasing imaging dose (# of hist.) from left to
right. Mass size of 5 mm and mass density of 1.1 remain constant.

Figure 6: Dose distribution (# of hist.) and percentages of optimal dose
considered by breast density.

Figure 7: Glandular dose distributions for the dataset.

Parameter | Considered Values
---|---
Breast phantom density | Dense, Hetero, Scattered, Fatty
Mass radius (mm) | 5, 7, 9
Mass density | 1.0, 1.06, 1.1
Relative Dose | 20$\%$, 40$\%$, 60$\%$, 80$\%$, 100$\%$
Detector type | DIR
Table 1: Parameters and their values corresponding to the M-SYNTH dataset.

##  4 Related Datasets

To date, a number of datasets for mammographic image analysis have been
collected (see Table 2). The majority of datasets are created from patient
data collected from DM [29, 30, 31, 32] or digital breast tomosynthesis (DBT)
[33, 34] scans from various clinical sites. The DREAM Challenge [35] offered
datasets for development of AI-based mammography analysis techniques. Patient
datasets vary widely in the types of labels available, and the data may be
biased toward the demographic characteristics of patients at the source site.
While there exist datasets, such as the EMory BrEast imaging Dataset (EMBED)
[34], that specifically focus on equal representation (in this case, equal
representation of African American and White patients), collecting a truly
balanced dataset across all possible characteristics may not be possible with
patient cases.

We found only two in silico datasets for mammography analysis. The first
dataset, published by Sarno [36], consists of 150 patient-derived digital
breast models with uncompressed computational breast phantoms derived from 3D
breast images acquired with an in-house dedicated breast computed tomography
(CT) scanner. The models were processed by a voxel classification algorithm
into four materials (air, adipose tissue, fibroglandular tissue, and skin).
The second dataset is the VICTRE [24] collection that consists of about 3,000
digital patients with breast sizes and densities representative of a screening
population. Digital microcalcification clusters and spiculated masses were
inserted in the voxelized phantoms to create the positive cohort. The phantoms
were imaged in silico to produce digital mammogram projections and digital
breast tomosynthesis volumes. In comparison to both of these datasets, our
work contains more significant variability in breast and mass characteristics,
as well as a range of applied dose levels for image acquisition, in order to
facilitate comparative evaluations of AI across characteristic changes.

Real patient datasets
---
Dataset | DM present | DBT present |  $\\#$ cases |  $\\#$ images | Image categories |  Population
Duke [33] | No | Yes | 5060 | 22032 | Cancer, benign, actionable, normal |  USA
ADMANI [29] | Yes | No | 629863 | 4411263a | Normal, recall |  Several countries
EMBED [34] | Yes | Yes | 116000 | 3383659b | Invasive cancer, non-invasive cancer, high risk, |  USAc
|  |  |  |  | borderline, benign, negative, non-breast cancer |
CMMD [30] | Yes | No | 1775 | 3712 | Benign, malignant |  China
INBreast [32] | Yes | No | 115 | 410 | Benign, malignant, normal |  Portugal
OPTIMAM [31] | Yes | No | 172,282 | 172,282 | Normal, interval cancers, benign, malignant |  UK
In silico datasets
Dataset | DM present | DBT present |  $\\#$ images | Image categories | Phantom variability |  Imaging
Sarno [36] | Yes | Yes | 150d | Normal | No |  No
VICTRE [24] | Yes | Yes | 2986 | Negative, positive cohort | Yese |  Yes
M-SYNTH (Ours) | Yes | Nog | 44914 | Negative, positive cohort | Yesf |  Yes

  * a

subset available for the RSNA Cancer Detection AI challenge.

  * b

20$\%$ available via AWS; contains annotated lesions.

  * c

equal representation of African american and White

  * d

150 uncompressed, 60 compressed images

  * e

four breast densities, same lesions across all positive cohort

  * f

3 lesion densities, 3 lesion sizes, 4 breast densities, 5 different doses

  * g

A corresponding DBT image dataset will be provided in a future release of the
dataset.

Table 2: Summary of existing mammographic image datasets.

##  5 Results and Analysis

In this section, we present an approach to using our M-SYNTH dataset to
evaluate an AI device. Formally, an image processing AI model $F$ takes as
input an image $r$ and predicts a specific property of interest $F(r)$ about
the image. For example, such a model can predict the presence or absence of a
mass. Typically for AI models, $F$ is a neural network and is trained on a
dataset of images and their labels
$T_{train}=\\{(r_{1},l_{1}),(r_{2},l_{2}),\ldots(r_{n},l_{n})\\}$, and then
evaluated on a held-out dataset $T_{test}$. When using patient images,
evaluation is limited to the variability contained in the samples and in the
annotations present across examples in the fixed test set $T_{test}$. Instead,
we propose to generate $T_{train}$ and $T_{test}$ dynamically using $D$ and
$I$ described above in order to test $F$ across variations in model $x$ and
acquisition parameters $y$.

###  5.1 Implementation Details

Evaluation Metrics We evaluate performance using the area under curve (AUC)
metric for a mass detection task. Specifically, we treat evaluation as a
multiple reader multiple case study, where an AI model is a single reader.
Multiple readers are obtained by re-training the model with different random
seeds. We rely on the iMRMC software [37, 38] to identify associated
confidence intervals.

Network Training We represent the AI-enabled device as a neural network with
an efficientnet_b0 architecture, receiving an image with one channel and
dimensions of 224 by 224, and outputting a binary mass presence label. The
network is trained with batch size 64 using binary cross entropy loss (BCE)
and optimized using RMSProp optimizer (with learning rate 0.0001). We rely on
the timm library [39] and fine-tune the model pre-trained with ImageNet [40].
We also compared performance with alternative architectures
(vit_small_patch16_224 and vgg_16), but results were very similar (see
supplementary material).

For each specific breast density, radiation dose level, and mass size and
density, the 300 images in the M-SYNTH dataset were divided into 200 for
training, 50 for validation, and 50 for test. For comparison, we also train
the AI device on 410 patient DM images from the INBreast dataset [32], where
images were obtained using MammoNovation Siemens full-field digital
mammography system with a solid-state amorphous selenium detector. We use the
same pre-processing and training regimes on this dataset and learn a network
to predict mass presence. The trained models on the real patient dataset were
then tested on 50 examples of M-SYNTH dataset for each specific breast
density, dose level, and mass size and density. The full experimental setup is
implemented in Python and C over a cluster with 50 Tesla V100-SXM2 GPUs.

###  5.2 Experimental Results

We identify two tasks that can be performed using our method. In the subgroup
analysis task, we train and test an AI model using the released synthetic
(M-SYNTH) dataset to identify performance changes on specified subgroups. In
the patient data evaluation task, we study how an AI model trained on patient
data (InBreast) performs on the proposed M-SYNTH dataset. This task can help
identify where the trained model may show variable performance for different
subgroups belonging to the target population.

(a) Mass size (b) Mass density (c) Breast density (d) Dose level

Figure 8: Subgroup analysis. Performance change across (a) mass size, (b) mass
density, (c) breast density, and (d) radiation dose, for models trained and
tested on our M-SYNTH dataset. These parameters remained constant for the set
of experiments performed during both training and test: (a) Fatty breast
phantom, mass density of 1.06, and relative dose of 100%. (b) Fatty breast
phantom, mass size of 7 mm, and relative dose of 100%. (c) Mass density of
1.06, mass size of 7 mm, and relative dose of 100%. (d) Fatty breast phantom,
mass density of 1.06, and mass size of 7 mm.

Subgroup Analysis. In Figures 8 and 9, we report the results of the AI model
performance at detecting masses, when the model is trained and tested on the
our dataset (see Section 5.1 for details of splits). We find that masses with
larger sizes or higher densities (Figures 8a-b) are more easily detected.
Although models trained on all sizes or mass densities have the highest
performance, when the models are trained on smaller masses or lower densities,
they generalize better to other masses (more difficult cases).The performance
of the models are highest when they are tested and trained on the same breast
density and decrease as the density of the test breast phantom differs from
the train phantom (Figures 8c). The dose levels applied in this study have
minimal impact on the performance of the models and resulted in similar AUC
values (Figures 8d). Evaluation of the performance change across all the
breast densities (Figures 9a-b) reveals that the AUC improves with larger mass
density and mass size, yet is impacted by the breast density, where mass
detection performance is lowest in high-density breasts (dense) and highest in
low-density breasts (fatty) in most of the cases, consistent with findings
from clinical practice.

(a) AUC as a function of mass size across all breast densities. (b) AUC as a
function of mass density across all breast densities.

Figure 9: Subgroup analysis. Performance changes for models trained and tested
on our M-SYNTH dataset. For each data point, the model is trained on 250
images with masses of radii of 7 mm and mass densities of 1.06, and tested on
50 images with mass characteristics shown in plots for each specific breast
density. The radiation dose level remains constant at 100% of the clinically
recommended dose for each breast density during training and test.

Patient Data Evaluation. In Figure 10, we report experiments where the AI
model is trained on INBreast data and evaluated on the M-SYNTH data. Although
the performance results for all experiments are lower in general, we find a
similar set of trends as when the model is trained on M-SYNTH data. Note that
we have made no attempt to match the radiation dose levels or the image
acquisition parameters for these comparisons using patient images. Even though
the simulated pipeline is designed to replicate a specific DM system with a
particular detector technology and technique factors, the comparison suggests
similarity between the datasets. The images are qualitatively different but
overall have similar glandular patterns which is an important consideration
for the realism of the task of detecting masses in a noisy background. We also
assessed similarity between INBreast and M-SYNTH datasets in terms of low-
level pixel distributions using first five statistical moments: mean,
variance, skewness, kurtosis, and hyperskewness. We found that there is a
reasonably good alignment in terms of moments, especially when the synthetic
images were included at all four breast densities (see supplementary
material). Future work should develop a more detailed comparison including
radiomics features for the training and testing datasets used in the study to
complement the validation of our approach.

(a) AUC as a function of mass size across all breast densities. (b) AUC as a
function of mass density across all breast densities.

Figure 10: Model Evaluation. Performance changes for a model trained on 410
real patient images (INBreast dataset) and tested on our M-SYNTH dataset. The
test sets consist of 50 images using parameters shown in the plots. The test
radiation dose is set to 100% of the clinically recommended dose for each
breast density.

Limitations. There are a number of limitations to our work. First, simulations
may require long runtimes and demand large computational resources, thus
somewhat limiting the amounts of data that can be generated. This limitation
needs to be considered with respect to the difficulty of obtaining large
patient image datasets with known mass locations. In addition, data can be
pre-generated offline (as we do with the M-SYNTH dataset), therefore, removing
the large runtime limit and computational burden off the user. Second, testing
with simulations is constrained to the variability captured by the parameter
space of the object models for anatomy and pathology and the acquisition
system. Thus, the complexity of the object model and acquisition system may
need to be adjusted depending on the complexity of the questions to be
investigated with simulated testing. In particular, a potential risk of
testing using simulated data is missing the variability observed in patient
populations. Finally, there is a risk of misjudging model performance due to a
domain gap between real and synthetic examples. However, the realism and
sophistication of object-based modeling of the imaging pipeline is improving
rapidly and may soon compete with other approaches, making approaches based on
synthetic data useful and practical for regulatory evaluation of AI-enabled
medical devices.

##  6 Conclusion and Future Work

We introduce and discuss an approach for validating AI models using physics-
based simulations of digital humans from the object space to the image data,
specifically for the task of breast cancer mass detection. The simulated
images are highly realistic and offer a challenging test case for AI model
evaluation. Our findings are consistent with expected performance and show
that the AI model performance increases with mass size and mass density as
expected. Finally, we show that our approach can be used to validate a model
trained on independent patient data. This finding suggests that the proposed
simulation setup can be used as a framework for more general evaluation of
medical AI devices. The goal of this study is to demonstrate as proof-of-
concept the feasibility of using simulated data to evaluate the comparative
performance of AI models. In future work, it would be important to assess the
evaluation approach for additional parameters in terms of the distribution of
the population of digital humans in the object space, and for a range of image
acquisition systems (e.g., by considering alternative simulators). By imaging
a more diverse population of breast models, we hope to identify additional
insights regarding AI evaluation. Finally, it is important to note that the
testing is limited to the variability captured in the digital representations
and may not fully indicate absolute real-world performance or trends. This
study illustrates that physics-based simulation of mammography images can
represent a less burdensome and cost-efficient approach for the evaluation of
AI model performance across a wide range of scenarios, including a variety of
image acquisition parameters and diverse populations that may not be available
or are hard to obtain from human studies. Moreover, this approach offers a
complementary evaluation paradigm that does not depend on the availability of
patient data.

##  7 Acknowledgements

We thank Andreu Badal (OSEL/CDRH/FDA) and anonymous reviewers for helpful
suggestions, Kenny Cha, Mike Mikailov and the OpenHPC team (OSEL/CDRH/FDA) for
providing help with experiments, Akhonda, Mohammad (OSEL/CDRH/FDA) for help
with data release, and Andrea Kim (OSEL/CDRH/FDA) for rendering visualizations
of the 3D breast model. This is a contribution of the US Food and Drug
Administration and is not subject to copyright. The mention of commercial
products herein is not to be construed as either an actual or implied
endorsement of such products by the Department of Health and Human Services.

## References

  * [1] A Badano, M Lago, E Sizikova, JG Delfino, S Guan, MA Anastasio, and B Sahiner.  The stochastic digital human is now enrolling for in silico imaging trials–methods and tools for generating digital cohorts.  arXiv preprint arXiv:2301.08719, 2023.
  * [2] Ana Barragán-Montero, Umair Javaid, Gilmer Valdés, Dan Nguyen, Paul Desbordes, Benoit Macq, Siri Willems, Liesbeth Vandewinckele, Mats Holmström, Fredrik Löfman, et al.  Artificial intelligence and machine learning for medical imaging: A technology review.  Physica Medica, 83:242–256, 2021.
  * [3] Christian G Graff.  A new, open-source, multi-modality digital breast phantom.  In Medical Imaging 2016: Physics of Medical Imaging, volume 9783, pages 72–81. SPIE, 2016.
  * [4] Andreu Badal, Diksha Sharma, Christian G Graff, Rongping Zeng, and Aldo Badano.  Mammography and breast tomosynthesis simulator for virtual clinical trials.  Computer Physics Communications, 261:107779, 2021.
  * [5] Ian Goodfellow, Jean Pouget-Abadie, Mehdi Mirza, Bing Xu, David Warde-Farley, Sherjil Ozair, Aaron Courville, and Yoshua Bengio.  Generative adversarial networks.  Communications of the ACM, 63(11):139–144, 2020.
  * [6] Youssef Skandarani, Pierre-Marc Jodoin, and Alain Lalande.  Gans for medical image synthesis: An empirical study.  arXiv preprint arXiv:2105.05318, 2021.
  * [7] Shuyue Guan and Murray Loew.  Breast cancer detection using synthetic mammograms from generative adversarial networks in convolutional neural networks.  Journal of Medical Imaging, 6(3):031411–031411, 2019.
  * [8] Yinhao Ren, Zhe Zhu, Yingzhou Li, Dehan Kong, Rui Hou, Lars J Grimm, Jeffery R Marks, and Joseph Y Lo.  Mask embedding for realistic high-resolution medical image synthesis.  In MICCAI, pages 422–430. Springer, 2019.
  * [9] Elena Sizikova, Xu Cao, Ashia Lewis, Kenny Moise, and Megan Coffee.  Improving computed tomography (ct) reconstruction via 3d shape induction.  arXiv preprint arXiv:2208.10937, 2022.
  * [10] Diederik P Kingma and Max Welling.  Auto-encoding variational bayes.  arXiv preprint arXiv:1312.6114, 2013.
  * [11] Danilo Rezende and Shakir Mohamed.  Variational inference with normalizing flows.  In International conference on machine learning, pages 1530–1538. PMLR, 2015.
  * [12] Piotr Bojanowski, Armand Joulin, David Lopez-Paz, and Arthur Szlam.  Optimizing the latent space of generative networks.  arXiv preprint arXiv:1707.05776, 2017.
  * [13] Jonathan Ho, Ajay Jain, and Pieter Abbeel.  Denoising diffusion probabilistic models.  Advances in Neural Information Processing Systems, 33:6840–6851, 2020.
  * [14] Florinel-Alin Croitoru, Vlad Hondru, Radu Tudor Ionescu, and Mubarak Shah.  Diffusion models in vision: A survey.  arXiv preprint arXiv:2209.04747, 2022.
  * [15] Weimin Zhou, Sayantan Bhadra, Frank J Brooks, Hua Li, and Mark A Anastasio.  Learning stochastic object models from medical imaging measurements by use of advanced ambient generative adversarial networks.  Journal of Medical Imaging, 9(1):015503, 2022.
  * [16] Jingkang Wang, Ava Pun, James Tu, Sivabalan Manivasagam, Abbas Sadat, Sergio Casas, Mengye Ren, and Raquel Urtasun.  Advsim: Generating safety-critical scenarios for self-driving vehicles.  In CVPR, pages 9909–9918, 2021.
  * [17] Naveed Akhtar and Ajmal Mian.  Threat of adversarial attacks on deep learning in computer vision: A survey.  Ieee Access, 6:14410–14430, 2018.
  * [18] Samuel G Finlayson, John D Bowers, Joichi Ito, Jonathan L Zittrain, Andrew L Beam, and Isaac S Kohane.  Adversarial attacks on medical machine learning.  Science, 363(6433):1287–1289, 2019.
  * [19] Hokuto Hirano, Akinori Minagi, and Kazuhiro Takemoto.  Universal adversarial attacks on deep neural networks for medical image classification.  BMC medical imaging, 21:1–13, 2021.
  * [20] Chaowei Xiao, Dawei Yang, Bo Li, Jia Deng, and Mingyan Liu.  Meshadv: Adversarial meshes for visual recognition.  In CVPR, pages 6898–6907, 2019.
  * [21] Xiaohui Zeng, Chenxi Liu, Yu-Siang Wang, Weichao Qiu, Lingxi Xie, Yu-Wing Tai, Chi-Keung Tang, and Alan L Yuille.  Adversarial attacks beyond the image space.  In CVPR, pages 4302–4311, 2019.
  * [22] Hsueh-Ti Derek Liu, Michael Tao, Chun-Liang Li, Derek Nowrouzezahrai, and Alec Jacobson.  Beyond pixel norm-balls: Parametric adversaries using an analytically differentiable renderer.  arXiv preprint arXiv:1808.02651, 2018.
  * [23] Guillaume Leclerc, Hadi Salman, Andrew Ilyas, Sai Vemprala, Logan Engstrom, Vibhav Vineet, Kai Xiao, Pengchuan Zhang, Shibani Santurkar, Greg Yang, et al.  3db: A framework for debugging computer vision models.  arXiv preprint arXiv:2106.03805, 2021.
  * [24] Aldo Badano, Christian G Graff, Andreu Badal, Diksha Sharma, Rongping Zeng, Frank W Samuelson, Stephen J Glick, and Kyle J Myers.  Evaluation of digital breast tomosynthesis as replacement of full-field digital mammography using an in silico imaging trial.  JAMA network open, 1(7):e185474–e185474, 11 2018.
  * [25] Steve A. Maas, Benjamin J. Ellis, Gerard A. Ateshian, and Jeffrey A. Weiss.  FEBio: Finite Elements for Biomechanics.  Journal of Biomechanical Engineering, 134(1):011005, 02 2012.
  * [26] Luis de Sisternes, Jovan G Brankov, Adam M Zysk, Robert A Schmidt, Robert M Nishikawa, and Miles N Wernick.  A computational model to generate simulated three-dimensional breast masses.  Medical physics, 42(2):1098–1118, 2015.
  * [27] Andrea Kim, Aunnasha Sengupta, and Aldo Badano.  Automated animation pipeline for visualizing in silico tumor growth models.  Physics of Medical Imging, 12463, 2023.
  * [28] Aunnasha Sengupta, Andreu Badal, Andrey Makeev, and Aldo Badano.  Computational models of direct and indirect x-ray breast imaging detectors for in silico trials.  Medical Physics, 49(11):6856–6870, 2022.
  * [29] Helen ML Frazer, Jennifer SN Tang, Michael S Elliott, Katrina M Kunicki, Brendan Hill, Ravishankar Karthik, Chun Fung Kwok, Carlos A Peña-Solorzano, Yuanhong Chen, Chong Wang, et al.  Admani: Annotated digital mammograms and associated non-image datasets.  Radiology: Artificial Intelligence, 5(2):e220072, 2022.
  * [30] Chunyan Cui, Li Li, Hongmin Cai, Zhihao Fan, Ling Zhang, Tingting Dan, Jiao Li, and Jinghua Wang.  The chinese mammography database (cmmd): An online mammography database with biopsy confirmed types for machine diagnosis of breast.  Data Cancer Imaging Arch, 2021.
  * [31] Mark D. Halling-Brown, Lucy M. Warren, Dominic Ward, Emma Lewis, Alistair Mackenzie, Matthew G. Wallis, Louise S. Wilkinson, Rosalind M. Given-Wilson, Rita McAvinchey, and Kenneth C. Young.  Optimam mammography image database: A large-scale resource of mammography images and clinical data.  Radiology: Artificial Intelligence, 3(1), 2020.
  * [32] Inês C Moreira, Igor Amaral, Inês Domingues, António Cardoso, Maria Joao Cardoso, and Jaime S Cardoso.  Inbreast: toward a full-field digital mammographic database.  Academic radiology, 19(2):236–248, 2012.
  * [33] Mateusz Buda, Ashirbani Saha, Ruth Walsh, Sujata Ghate, Nianyi Li, Albert Święcicki, Joseph Y Lo, and Maciej A Mazurowski.  Dete ction of masses and architectural distortions in digital breast tomosynthesis: a publicly available dataset of 5,060 patients and a deep learning model.  arXiv preprint arXiv:2011.07995, 2020.
  * [34] Jiwoong J Jeong, Brianna L Vey, Ananth Bhimireddy, Thomas Kim, Thiago Santos, Ramon Correa, Raman Dutt, Marina Mosunjac, Gabriela Oprea-Ilies, Geoffrey Smith, et al.  The emory breast imaging dataset (embed): A racially diverse, granular dataset of 3.4 million screening and diagnostic mammographic images.  Radiology: Artificial Intelligence, 5(1):e220047, 2023.
  * [35] Berkman Sahiner.  Digital mammography dream challenge overview (conference presentation).  In Medical Imaging 2017: Computer-Aided Diagnosis, volume 10134, pages 1159–1159. SPIE, 2017.
  * [36] Antonio Sarno, Giovanni Mettivier, Francesca di Franco, Antonio Varallo, Kristina Bliznakova, Andrew M. Hernandez, John M. Boone, and Paolo Russo.  Dataset of patient-derived digital breast phantoms for in silico studies in breast computed tomography, digital breast tomosynthesis, and digital mammography.  Medical Physics, 48(5):2682–2693, 2021.
  * [37] Brandon D. Gallas, Andriy Bandos, Frank W. Samuelson, and Robert F. Wagner.  A framework for random-effects roc analysis: Biases with the bootstrap and other variance estimators.  Communications in Statistics - Theory and Methods, 38(15):2586–2603, 2009.
  * [38] RST Catalog.  iMRMC: Software for the Statistical Analysis of multi-reader multi-case studies, June 2022.
  * [39] Ross Wightman.  Pytorch image models.  https://github.com/rwightman/pytorch-image-models, 2019.
  * [40] Jia Deng, Wei Dong, Richard Socher, Li-Jia Li, Kai Li, and Li Fei-Fei.  Imagenet: A large-scale hierarchical image database.  In CVPR, pages 248–255. Ieee, 2009.

##  8 Supplementary Material

###  8.1 Data Availability

M-SYNTH and code for processing can be found in
https://github.com/DIDSR/msynth-release. Please following the instructions on
Github to dowload files from Huggingface. M-SYNTH is organized into a
directory structure that indicates the parameters. The folder

    data/device_data_VICTREPhantoms_spic_[LESION_DENSITY]/[DOSE]/[BREAST_DENSITY]/
    2/[LESION_SIZE]/SIM/P2_[LESION_SIZE]_[BREAST_DENSITY].8337609.[PHANTOM_FILE_ID]/
    [PHANTOM_FILEID]/

contains image files imaged with the specified parameters. Each folder
contains mammogram data that can be read from .raw format (.mhd contains
supporting data), or DICOM (.dcm) format. Note that only examples with odd
`PHANTOM_FILEID` contain lesions, others do not. Coordinates of lesions can be
found in .loc files. For instance:

    --P2_5.0_hetero.8337609.1/1/
    ----DICOM_dm
    ------000.dcm
    ----projection_DM1.loc
    ----projection_DM1.mhd
    ----projection_DM1.raw

contains a lesion-present breast example with mass size (radius) of 5.0 mm
(approximate, as the mass is not perfectly spherical), mass density 1.0, dose
(# histories) $1.02\times 10^{10}$, and heterogeneously dense breast density.
Code and dataset is released with the Creative Commons 1.0 Universal License
(CC0).

###  8.2 Timing Analysis

We now review the timing required to perform mass insertion and imaging.
Timings were computed on a Tesla V100-PCIE GPU card with 32 GB RAM. In Table
3, we review the mean timing (in minutes) for mass insertion by breast density
and mass size across each category of examples. We find that larger mass size
requires a slight increase in time. However, breast density significantly
affects timing because the reading and writing times are proportional to the
number of voxels in the volume. In particular, lower density breasts, which
are larger in size on the average, need more insertion time, with fatty
breasts requiring nearly 3.5 as much time than dense breasts. Note that mass
density is set during projection, therefore, it does not affect insertion
time.

Breast Density | Mass Size (mm) | Time (min)
---|---|---
Fatty | 5.0 | 7.152661
| 7.0 | 7.206867
| 9.0 | 7.337922
Scattered | 5.0 | 5.035144
| 7.0 | 5.139315
| 9.0 | 5.366446
Hetero | 5.0 | 2.583082
| 7.0 | 2.769962
Dense | 5.0 | 2.095512
| 7.0 | 2.327806
Table 3: Timing analysis for mass insertion by breast density and mass size.

In Table 4, we review the imaging time required for each breast density. The
time varies from 2.84 min for most dense to 13.46 min to least dense breasts.
Note that total time for creating of each DM image is either the imaging time
(no mass inserted) or imaging + mass insertion times. Given our high
performance cluster with access to multiple GPUs (where each example requires
access to one GPU), we were able to generate the complete dataset in about two
weeks.

Breast Density | Time (min)
---|---
Fatty | 13.463809
Scattered | 11.002291
Hetero | 3.655613
Dense | 2.842028
Table 4: Timing analysis for imaging by breast density.

###  8.3 Rendering of Breast Phantoms

Additional renderings of the breast phantoms generated for the study are shown
in Figure 11, demonstrating a high level of detail and anatomical variability
within and among models.

Figure 11: Renderings of the breast phantoms for each composition.

###  8.4 Real and Synthetic Image Similarity Assessment

In order to investigate the similarity in terms of low-level pixel
distributions between the real patient (INBreast) and synthetic (M-SYNTH)
datasets, we estimated the first five statistical moments (mean, variance,
skewness, kurtosis, and hyperskewness). Although there is a differences
between synthetic and real examples, the distributions and ranges are
reasonably aligned.

Figure 12: First five statistical moments for the real patient (INBreast, 410
images) and synthetic (M-SYNTH, 1200 images consisted of 300 images for each
breast density) datasets. The measurements were performed on images with mass
size of 7 mm, mass density of 1.06, and at 100% clinically recommended dose.
m0: mean, m1: variance, m2: skewness, m3: kurtosis, and m4: hyperskewness.

###  8.5 Additional Subgroup Analysis

####  8.5.1 Mass Size and Density Effects

We further study the impact of generalization of the training dataset on the
performance of mass detection. In Figure 13a, we train the models on
individual mass sizes, as well as on all the sizes. The training mass density
of 1.06 and relative radiation dose of 100% are kept constant. Each model is
trained and tested on the same breast density that is given on top of each
figure, with the test mass size and mass density as shown. We find that the
models trained on all sizes (dashed lines) have an equal or better performance
on small masses (i.e., 5 mm) than the models trained on a specific mass radius
(solid lines) (except for scattered breast density). However the models
trained on all sizes generalize worse to the larger masses, compared to the
models trained and tested on the same mass size. Similarly, in Figure 13b, we
train the models on individual mass densities, as well as on all the mass
densities. The training mass size of 7 mm and relative radiation dose of 100%
are kept constant. Each model is trained and tested on the same breast density
that is given on top of each figure, with the test mass density and mass size
as shown. We find that in most of the cases, the models trained on all the
mass densities (dashed lines) result in worse performance than the models
trained on a specific mass density (solid lines), specially as the test mass
size increases. Thus, these models are not able to generalize well to masses
with different densities on the testing dataset.

(a) AUC as a function of mass size across all breast densities. (b) AUC as a
function of mass density across all breast densities.

Figure 13: Performance changes for models trained and tested on our M-SYNTH
dataset. For each data point, the model is trained on 250 images with (a)
masses of radii of 7 mm and mass densities of 1.06 (solid lines, —) or all
mass densities (dashed lines, - - - ), (b) masses of radii of 7 mm (solid
lines, —) or all sizes (dashed lines, - - - ) and mass densities of 1.06. The
model is tested on 50 images with mass characteristics shown in plots for each
specific breast density. The radiation dose level remains constant at 100% of
the clinically recommended dose for each breast density during training and
test.

####  8.5.2 Network Architecture Effects

In order to evaluate the effect of the AI enabled device, we repeat the
experiments with additional model architectures of vit_small_patch16_224 and
vgg_16. As shown in Figures 14 and 15, using different models results in
similar results and has minimal impact of the outcome of the experiments.

(a) AI enabled device architecture: efficientnet_b0 (b) AI enabled device
architecture: vit_small_patch16_224 (c) AI enabled device architecture: vgg_16

Figure 14: Performance changes as a function of mass size across all breast
densities. Different architectures of (a) efficientnet_b0, (b)
vit_small_patch16_224, and (c) vgg_16 are used as the AI enabled device to be
trained and tested on our M-SYNTH dataset. For each data point, the model is
trained on 250 images with masses of radii of 7 mm and mass densities of 1.06,
and tested on 50 images with mass characteristics shown in plots for each
specific breast density. The radiation dose level remains constant at 100% of
the clinically recommended dose for each breast density during training and
test.

(a) AI enabled device architecture: efficientnet_b0 (b) AI enabled device
architecture: vit_small_patch16_224 (c) AI enabled device architecture: vgg_16

Figure 15: Performance changes as a function of mass density across all breast
densities. Different architectures of (a) efficientnet_b0, (b)
vit_small_patch16_224, and (c) vgg_16 are used as the AI enabled device to be
trained and tested on our M-SYNTH dataset. For each data point, the model is
trained on 250 images with masses of radii of 7 mm and mass densities of 1.06,
and tested on 50 images with mass characteristics shown in plots for each
specific breast density. The radiation dose level remains constant at 100% of
the clinically recommended dose for each breast density during training and
test.


# Discovery of a Neutron Star Oscillation Mode During a Superburst

Tod Strohmayer1 and Simin Mahmoodifar2
1Astrophysics Science Division and Joint Space-Science Institute, NASA’s
Goddard Space Flight Center, Greenbelt, MD 20771, USA
2Department of Physics and Joint Space-Science Institute, University of
Maryland College Park, MD 20742, USA

###### Abstract

Neutron stars are among the most compact objects in the universe and provide a
unique laboratory for the study of cold ultra-dense matter. While
asteroseismology can provide a powerful probe of the interiors of stars, for
example, helioseismology has provided unprecedented insights about the
interior of the sun, comparable capabilities for neutron star seismology have
not yet been achieved. Here we report the discovery of a coherent X-ray
modulation from the neutron star 4U 1636$-$536 during the February 22, 2001
thermonuclear superburst seen with NASA’s Rossi X-ray Timing Explorer (RXTE)
that is very likely produced by a global oscillation mode. The observed
frequency is $835.6440\pm 0.0002$ Hz (1.43546 times the stellar spin frequency
of 582.14323 Hz) and the modulation is well described by a sinusoid
($A+B\sin(\phi-\phi_{0})$) with fractional half-amplitude of $B/A=0.19\pm
0.04\%$ ($4-15$ keV). The observed frequency is consistent with the expected
inertial frame frequency of a rotationally-modified surface g-mode, an
interfacial mode in the ocean-crust interface or perhaps an r-mode. Observing
an inertial frame frequency–as opposed to a co-rotating frame
frequency–appears consistent with the superburst’s thermal emission arising
from the entire surface of the neutron star, and the mode may become visible
by perturbing the local surface temperature. We briefly discuss the
implications of the mode detection for the neutron star’s projected velocity
and mass. Our results provide further strong evidence that global oscillation
modes can produce observable modulations in the X-ray flux from neutron stars.

###### Subject headings:

stars: neutron — stars: oscillations — stars: rotation — X-rays: binaries —
X-rays: individual (4U 1636$-$536) — methods: data analysis

##  1\. Introduction

Neutron stars provide natural laboratories for the study of ultra-dense
matter. A primary method for such studies is to accurately measure the mass-
radius relationship for neutron stars, which depends directly on the equation
of state of dense matter. However, different phases of ultra-dense matter can
have similar equations of state, such that mass-radius measurements alone may
not be definitive in probing the composition of matter at the highest
densities (Alford et al., 2013; Lattimer, 2012). Additional, complementary
observables sensitive to the phase of dense matter would be extremely valuable
for a more comprehensive understanding of neutron star interiors. The
frequencies and damping timescales of different global oscillation modes of
neutron stars, such as the Rossby waves (r-modes) and gravity modes (g-modes),
which have the Coriolis force and buoyancy as their respective restoring
forces, depend on their interior structure and composition (Yoshida & Lee,
2000; Alford et al., 2012). Thus, observations of the oscillation modes of a
neutron star can, in principle, provide such a complementary probe of its
interior properties (Passamonti & Andersson, 2011).

The low mass X-ray binary (LMXB) 4U 1636$-$536 (hereafter, 4U1636) is a well
studied accreting neutron star binary with an optically-determined orbital
period of 3.8 hr, and is also a well-known X-ray burst source (Augusteijn et
al., 1998; Giles et al., 2002). Oscillations at $\approx 582$ Hz are detected
during some of its thermonuclear bursts and are almost certainly spin
modulation pulsations (Strohmayer et al., 1998). This conclusion is further
supported by the detection of $582$ Hz pulsations during a portion of the
thermonuclear superburst observed from this source on 22 February, 2001 with
the Proportional Counter Array (PCA) onboard RXTE (Strohmayer & Markwardt,
2002). Superbursts are rare, energetic X-ray flares observed from LMXBs that
are likely caused by unstable thermonuclear burning of a carbon-rich layer
formed from the ashes of normal Type I X-ray bursts. The time evolution of the
582 Hz pulsations during the superburst was entirely consistent with the
system’s known orbital ephemeris, and modeling of the frequency drift enabled
some constraints on the neutron stars’s projected orbital velocity and its
rest-frame spin frequency (Strohmayer & Markwardt, 2002; Casares et al.,
2006).

Figure 1.—  Light curve of the portion of the superburst from 4U1636 used in
our pulsation search, along with the frequency evolution from our best-fitting
orbit model (red curve, right axis). The data are the $4-15$ keV PCA count
rates (left axis) in 1 s bins. The gap in the range from $\approx 2500-4000$ s
results from Earth occultation of the RXTE satellite. The first data interval
begins just prior to the peak of the superburst. Time zero is MJD
51962.70863263 (referred to the TDB timescale).

The bright, hours-long superburst from 4U1636 provides a unique opportunity to
search for X-ray modulations produced by global neutron star oscillation modes
for several reasons. First, the combination of its high luminosity (and thus
high counting rate in the PCA) and long duration provides excellent
sensitivity to potentially weak modulations. Second, the neutron star orbit
constraint derived from the 582 Hz pulsations enables a coherent search by
removal of the phase delays due to the neutron star’s orbital motion
(Strohmayer & Mahmoodifar, 2014). Finally, the thermonuclear burning provides
a plausible mechanism for the excitation of modes (the so-called
$\epsilon$-mechanism), and the shock wave produced by detonation of the carbon
fuel which powers the superburst is also likely to excite wave motions at or
near the surface (McDermott & Taam, 1987; Strohmayer & Lee, 1996; Piro &
Bildsten, 2004; Keek et al., 2012).

Figure 2.—  A portion of the power spectrum (Leahy-normalized) computed from
the 8000 s superburst interval shown in Figure 1. The x-axis shows Fourier
frequency scaled by the pulsar spin frequency. The plotted range includes both
the co-rotating ($0.4166\leq\sigma/\Omega\leq 0.7567$) and inertial
($1.243\leq\sigma/\Omega\leq 1.583$) frame frequency ranges searched. The
pulsar signal appears at 1 in these units (labeled), and the oscillation mode
peak is evident near 1.43546, with a Leahy power close to 50. The dashed
horizontal line denotes the $3\sigma$ detection level assuming a number of
trials equal to $4\times 3.17\times 10^{6}$. The factor of 4 conservatively
accounts for the number of energy ranges searched.

##  2\. A Coherent Mode Search During the Superburst

We carried out a coherent pulsation search during the superburst in the manner
described in detail in Strohmayer & Mahmoodifar (2014). As the methods are
clearly described there we omit some of the details here. We used the high
time resolution event mode (E_125us_64M_0_1s) data obtained for most of the
first two RXTE orbits in which the superburst was observed (see Strohmayer &
Markwardt 2002 for a discussion of the PCA data obtained during the
superburst). For the orbit correction we began with the values that minimized
a $\chi^{2}$ fit to the phase timing residuals of the 582 Hz pulsation, but
with the orbit period fixed at the well-known value of 0.15804693(16) days
(Strohmayer & Markwardt 2002; Giles et al. 2002). This gave a pulsar frequency
of $582.14323$ Hz, a projected neutron star velocity of 134.53 km/s, and an
epoch of $T_{90}=51962.74302977$ MJD (referred to Barycentric Dynamical Time,
TDB). These values are all consistent with those reported by Strohmayer &
Markwardt (2002). After correcting the event times for the neutron star’s
orbital motion, we created light curves sampled at 4096 Hz spanning 8000 s of
the portions of the first two RXTE orbits with high time resolution data (see
Figure 1 in Strohmayer & Markwardt 2002). We generated light curves in four
energy bands, beginning with the full PCA band, and then restricted the range
to more closely match that of the superburst’s thermal emission component
(Keek et al., 2014). Our last and tightest energy range spans $\approx 4-15$
keV, where most of the detected superburst thermal emission falls. This light
curve re-sampled to 1 second bins is shown in Figure 1 along with the orbital
frequency model used to correct the event times (red curve).

We then computed power spectra and searched the same frequency ranges as
described in Strohmayer & Mahmoodifar (2014). These encompass the expected
range of frequencies–relative to the stellar spin frequency–for r-modes and
some rotationally modified g-modes, if a mode modulates neutron star emission
at either the mode’s co-rotating ($0.4166<\sigma/\Omega<0.7567$) or inertial
frame frequency ($1.243<\sigma/\Omega<1.583$). Here, $\sigma$ is a Fourier
frequency and $\Omega$ is the stellar spin frequency. We note that for the
superburst’s thermal emission component–which almost certainly reflects
emission from the whole surface of the neutron star–it would seem the inertial
frame frequencies are the more likely to be observed (Heyl, 2005; Lee &
Strohmayer, 2005). This is in contrast to the case of thermal emission from
the hot-spot of an accreting millisecond X-ray pulsar (AMXP), where the mode
can periodically perturb the shape of the hot-spot that is fixed in the
rotating frame, producing an observed modulation at the mode’s co-rotating
frame frequency (Numata & Lee, 2010).

Our analysis revealed a significant power spectral peak at a frequency of
$835.644$ Hz ($\sigma/\Omega=1.43546$), in the inertial frame frequency search
range. The peak has a Leahy-normalized power of 49.3, and we estimate its
significance as $\exp(-49.3/2)\times N_{trials}$, where $N_{trials}=4\times
3.17\times 10^{6}$ is the total number of Fourier frequency bins in the two
ranges searched. The factor of 4 conservatively accounts for the number of
energy channel ranges we used to compute light curves and power spectra, and
results in a significance of $2.5\times 10^{-4}$. The estimate is conservative
because the channel ranges searched were not all independent. We also carried
out extensive Monte Carlo simulations to assess the significance, taking into
account the non-independence of the energy channel ranges searched. From these
simulations we find a significance of $\approx 1.5\times 10^{-4}$, which gives
us additional confidence in the detection. Figure 2 shows a portion of the
full power spectrum ($0.4\leq\sigma/\Omega\leq 1.6$) which includes the two
frequency ranges searched, as well as the range around the known pulsar
frequency. There are clearly two significant peaks present in the power
spectrum, one is the known pulsar signal (at 1 in these units), and the other
is the putative oscillation mode frequency at $1.43546$. There are no other
significant peaks detected. As a further check, we used the Fourier frequency
bins above $\sigma/\Omega=1.6$ to investigate the noise power distribution,
and find that it closely matches the expected $\exp(-P/2)$ distribution,
giving us further confidence that our significance estimates are robust.

The 835 Hz signal we detect is consistent with a coherent oscillation over the
8000 s of the light curve. We emphasize that prior searches, such as those
which discovered the 582 Hz pulsations, would not have detected this
modulation because they did not account for the motion of the neutron star,
which smears a narrow band signal into many adjacent Fourier frequency bins.
An estimate of this effect can be gleaned from Figure 1. The full (peak to
trough) fractional change in the frequency due to the orbital motion is
$\Delta\nu/\nu\approx 0.5/582=8.6\times 10^{-4}$. A coherent (single bin)
signal at 835 Hz will be smeared into $N_{f}=835\times(\Delta\nu/\nu)/\delta
f$ bins, where $\delta f$ is just the frequency resolution of the power
spectrum, in this case $\delta f=1/8000\;{\rm s}$, and
$N_{f}=835*(0.5/582)*8000\approx 5700$. Any such signal would be swamped by
the poisson noise in this large number of Fourier bins (Wood et al., 1991).

Figure 3.—  Pulse profile in the $\approx 4-15$ keV band obtained by folding
the light curve data at the oscillation mode frequency of
$\sigma/\Omega=1.43546$ in 12 phase bins. Two cycles are plotted for clarity.
The best fitting sinusoidal model (solid curve), $A+B\sin(\phi-\phi_{0})$, is
also plotted. The fractional half-amplitude of the best fitting model is
$B/A=0.19\pm 0.04\%$.

We folded the light curve data ($4-15$ keV) at the detected frequency and fit
a sinusoidal model of the form $A+B\sin(\phi-\phi_{0})$ to the resulting
profile (see Figure 3). This model provides an excellent fit to the data, and
we find a fractional half-amplitude of $B/A=0.19\pm 0.04\%$. We do not find
any evidence for significant harmonics associated with this signal. Assuming
the observed frequency is the inertial frame frequency of an oscillation mode
($\omega_{i}/\Omega=1.43546$), we can convert it to a co-rotating frame
frequency if we know its azimuthal wavenumber, $m$, since the two frequencies
are related as $\omega_{i}=m\Omega-\omega$. Interestingly, for $m=2$ we find
that $\omega/\Omega=2-1.43546=0.56454$, which is close to the candidate
frequency ratio of 0.57276 we identified in the AMXP XTE J1751$-$305
(Strohmayer & Mahmoodifar 2014), suggesting that the frequencies identified in
the two sources could perhaps be associated with the same oscillation mode.
For the AMXP source the modulation mechanism is more likely due to
oscillation-induced perturbations to the hot-spot fixed in the rotating frame,
whereas for the superburst, emission from the entire surface can be modulated
by a mode, perhaps due to local variations in the surface temperature. As
mentioned above, these processes would naturally lead to modulations at the
co-rotating and inertial frame frequencies, respectively.

##  3\. Discussion

Recent theoretical work has shown that the mode frequency ratio relative to
the spin frequency indicated for J1751$-$305 is consistent, in principle, with
the $l=m=2$ r-mode for neutron stars with plausible masses and radii
(Andersson et al., 2014), and surface g-modes with $l=1$ or 2 (Strohmayer &
Mahmoodifar 2014). To the extent that an $m=2$ r-mode identification is
appropriate for the modulation during the superburst, and since the mode
frequency ratios are similar, the same r-mode identification appears possible
for 4U1636 as well. The $l=m=2$ r-mode is the most unstable to the emission of
gravitational waves and its amplitude can grow exponentially–potentially
causing a rapid spin-down of the star–if viscosity and other non-linear
dissipation processes are not strong enough (Andersson, 1998; Lindblom et al.,
1998; Owen et al., 1998). The amplitude of the r-modes determines how fast
they can spin down the star. Now, the r-mode amplitude likely required to
account for the observed X-ray modulation during the superburst would, if
continuously present, produce a strong spin-down of the star (Andersson et al.
2014; Strohmayer & Mahmoodifar 2014). Such a spin-down is inconsistent with
the pulse timing data for J1751$-$305 (Strohmayer & Mahmoodifar 2014;
Andersson et al. 2014), but we note that there are no comparable long term
spin-down measurements for 4U1636, as the source has never been detected as a
persistent pulsar. Excitation of the mode to an observable amplitude might
also be intermittent, perhaps driven by the superburst’s energy release.

However, these inferences do not take account of the presence of a solid crust
within the neutron star. It has recently been argued that the r-mode amplitude
may be “amplified” in the surface layers above the crust compared to the core
because the crust can isolate the surface layers from the core motions (Lee,
2014). Indeed, recent estimates of r-mode amplitudes–at least those in the
core–suggest a low saturation amplitude of $10^{-7}-10^{-8}$ (Haskell et al.,
2012; Mahmoodifar & Strohmayer, 2013). Amplification of the mode in the
surface layers might then account for the observed X-ray modulations, and if
some process limits the r-mode saturation amplitude in the core (Bondarescu &
Wasserman, 2013), then the apparent inconsistency with the long-term spin
evolution in J1751$-$305 could perhaps be mitigated. Alternatively, surface
g-mode or ocean-crust interface mode identifications remain viable, and we
note that rotational “squeezing” of the mode displacements toward the
rotational equator does not necessarily pose a problem for the visibility of
the mode in either of these cases, as the superburst involves emission from
the entire neutron star surface, and not from a hot-spot localized near the
rotational pole (Bildsten et al., 1996; Piro & Bildsten, 2005).

Figure 4.—  Constraints on the component masses of 4U1636 assuming the neutron
star velocity of $134.5\pm 0.2$ km/s suggested by our mode frequency
detection. The orange region is excluded as zero age main sequence (ZAMS)
stars above this limit will have radii larger than the Roche lobe radius for
4U1636. The blue region below the “eclipse line” is excluded as these systems
would not be able to satisfy the velocity constraint (taken as $v_{ns}\geq
134.3$ km/s) at an inclination angle which would preclude eclipses. The
maximum allowed mass under these assumptions is $M_{ns}\leq 1.108M_{\odot}$,
and is marked by the vertical dashed line. The region of allowed masses for a
system inclination angle of $i=65^{\circ}$ (inside the parallel blue lines) is
also shown, to indicate the precision which could be attained if both the
inclination and velocity were accurately known.

If the putative inertial frame frequency at $\omega_{i}/\Omega=1.43546$ is
indeed constant over the 8000 s light curve interval, then it suggests that
the neutron star velocity can be much more tightly constrained than indicated
by the time evolution of the 582 Hz pulsations alone (Strohmayer & Markwardt
2002; Casares et al. 2006). That is because the 582 Hz pulsations were
detected over a much shorter fraction of the orbital period (800 s) than the
8000 s detection interval of the 835 Hz modulation. The much shorter pulse
train of the 582 Hz oscillation was consistent with a much larger range of
neutron star orbital velocity ($90-175$ km/s) than the 835 Hz signal
(Strohmayer & Markwardt 2002). Indeed, varying the orbital parameters away
from the values for which the $835.644$ Hz signal was detected results in a
“smearing-out” of the signal peak. As noted toward the end of §2, this is the
expected behavior of a coherent modulation when the chosen orbital parameters
do not accurately match the “true” values.

Assuming the mode frequency is constant in the star’s frame, and using it as a
“matched filter,” we can, in principle, constrain the projected neutron star
velocity to $134.5\pm 0.2$ km/s. A velocity outside of this range results in
significant signal loss in the $835.644$ Hz oscillation. We caution, however,
that this is a somewhat indirect means to infer the star’s velocity because we
do not have independent verification that the mode frequency is constant in
the star’s frame. For example, some intrinsic variation of the mode frequency
could mean that the orbital velocity which maximizes the 835 Hz signal is
different from the true orbital velocity. Previous radial velocity studies of
4U1636 based on both optical spectroscopy and the 582 Hz superburst pulsations
suggested a lower projected neutron star velocity, in the range from $90-113$
km/s. This was based on a “spectroscopic ephemeris” for the epoch (phase) of
superior conjunction of the companion deduced from Doppler tomography of the
Bowen emission lines observed from the optical companion (V801 Ara; Casares et
al. 2006). The offset in the “spectroscopic” and “photometric” phases (used
here in our $835.643$ Hz signal detection) corresponds to $\approx 2.5\sigma$
of the statistical uncertainty given for the reference epoch of the
spectroscopic ephemeris (Casares et al. 2006), however, the inclusion of
systematic uncertainties associated with identifying accurately the site of
Bowen fluorescence in the binary, as well as the (smaller) uncertainty in the
orbital period, would further reduce this potential discrepancy. We also note
that the timing analyses of the 582 Hz pulsations modestly favor—in yielding
lower $\chi^{2}$ timing residuals—the higher neutron star velocity suggested
by our $835.643$ Hz oscillation detection (Strohmayer & Markwardt 2002).

Nevertheless, such a high neutron star velocity would have important
implications for its mass, as we now explain. Constraints on the component
masses of 4U1636 assuming a neutron star velocity of $134.5\pm 0.2$ km/s are
shown in Figure 4. Systems within the orange shaded region are disfavored as
any zero age main sequence star (ZAMS) is larger than the Roche lobe radius
for 4U1636. We used a standard estimate for the radius of the Roche lobe
(Eggleton, 1983) and the stellar mass-radius relation of Tout et al. (1996) to
determine the lower boundary of this region. The “eclipse line” denotes
systems which give a neutron star velocity equal to the $1\sigma$ lower bound
of $134.3$ km/s at an inclination angle for which eclipses just begin. Since
eclipses are not seen from 4U1636, the parameter space below this line is
excluded. We used Kepler’s law and solved the Roche lobe eclipse geometry for
a point source (Chanan et al., 1976) to draw this line. The allowed range of
masses is also shown if the system inclination were $i=65^{\circ}$ (within the
parallel blue lines), and is meant to suggest the accuracy that could be
achieved if both the velocity and inclination were well constrained. The
vertical dotted line marks the maximum neutron star mass under these
assumptions, $M_{ns}\leq 1.108M_{\odot}$. From Figure 4 it is evident that an
orbital velocity as high as $134.5\pm 0.2$ km/s would strongly favor a low
mass neutron star in 4U1636. Interestingly, there is growing evidence for
“light” neutron stars, including the X-ray binaries 4U 1538$-$52 and SMC X-1
(Rawls et al., 2011; Coe et al., 2013). As discussed by Lattimer (2012), the
physics of core-collapse and the thermodynamics of the resulting lepton-rich
proto-neutron star, suggests a minimum mass of about $0.9-1.1$ $M_{\odot}$
depending on the entropy profile of the star (Strobel et al., 1999). Thus,
accurate measurements near the minimum would provide important insights into
the core-collapse physics (Lattimer, 2012). Finally, we emphasize that an
independent confirmation of the neutron star’s orbital velocity is crucial to
definitively establish such a low mass for 4U1636.

##  4\. Summary and Conclusions

We carried out a search for high frequency coherent X-ray modulations during
the February 2001 superburst from the LMXB 4U 1636$-$536 observed with RXTE.
We used a circular orbit model deduced from the time evolution of the 582 Hz
pulsations seen previously during the superburst to remove the Doppler delays
due to binary motion of the neutron star. We detected a coherent modulation at
a frequency of $835.6440\pm 0.0002$ Hz (1.43546 times the stellar spin
frequency of 582.14323 Hz), with a significance estimated from Monte Carlo
simulations of $\approx 1.5\times 10^{-4}$, and a fractional half-amplitude of
$0.19\pm 0.04\%$. This frequency is consistent, in principle, with the
expected inertial frame frequency of an $l=m=2$ r-mode, a rotationally
modified surface g-mode, or an ocean-crust interfacial mode. Our results
provide further strong evidence that global oscillation modes can produce
observable modulations in the X-ray flux from neutron stars. Sensitive
searches for such signals should be a priority for future, large area X-ray
timing missions.

We thank Cole Miller, Tony Piro and the anonymous referee for helpful comments
and discussions. TS acknowledges NASA’s support for high energy astrophysics.
SM acknowledges the support of the U.S. Department of Energy through grant
number DEFG02- 93ER-40762.

## References

  * Alford et al. (2013) Alford, M. G., Han, S., & Prakash, M. 2013, Phys. Rev. D, 88, 083013
  * Alford et al. (2012) Alford, M. G., Mahmoodifar, S., & Schwenzer, K. 2012, Phys. Rev. D, 85, 024007
  * Andersson (1998) Andersson, N. 1998, ApJ, 502, 708
  * Andersson et al. (2014) Andersson, N., Jones, D. I., & Ho, W. C. G. 2014, MNRAS, 442, 1786
  * Augusteijn et al. (1998) Augusteijn, T., van der Hooft, F., de Jong, J. A., van Kerkwijk, M. H., & van Paradijs, J. 1998, A&A, 332, 561
  * Bildsten et al. (1996) Bildsten, L., Ushomirsky, G., & Cutler, C. 1996, ApJ, 460, 827
  * Bondarescu & Wasserman (2013) Bondarescu, R., & Wasserman, I. 2013, ApJ, 778, 9
  * Casares et al. (2006) Casares, J., Cornelisse, R., Steeghs, D., et al. 2006, MNRAS, 373, 1235
  * Chanan et al. (1976) Chanan, G. A., Middleditch, J., & Nelson, J. E. 1976, ApJ, 208, 512
  * Coe et al. (2013) Coe, M. J., Angus, R., Orosz, J. A., & Udalski, A. 2013, MNRAS, 433, 746
  * Eggleton (1983) Eggleton, P. P. 1983, ApJ, 268, 368
  * Giles et al. (2002) Giles, A. B., Hill, K. M., Strohmayer, T. E., & Cummings, N. 2002, ApJ, 568, 279
  * Haskell et al. (2012) Haskell, B., Degenaar, N., & Ho, W. C. G. 2012, MNRAS, 424, 93
  * Heyl (2005) Heyl, J. S. 2005, MNRAS, 361, 504
  * Keek et al. (2014) Keek, L., Ballantyne, D. R., Kuulkers, E., & Strohmayer, T. E. 2014, ApJ, 789, 121
  * Keek et al. (2012) Keek, L., Heger, A., & in’t Zand, J. J. M. 2012, ApJ, 752, 150
  * Lattimer (2012) Lattimer, J. M. 2012, Annual Review of Nuclear and Particle Science, 62, 485
  * Lee (2014) Lee, U. 2014, MNRAS, 442, 3037
  * Lee & Strohmayer (2005) Lee, U., & Strohmayer, T. E. 2005, MNRAS, 361, 659
  * Lindblom et al. (1998) Lindblom, L., Owen, B. J., & Morsink, S. M. 1998, Physical Review Letters, 80, 4843
  * Mahmoodifar & Strohmayer (2013) Mahmoodifar, S., & Strohmayer, T. 2013, ApJ, 773, 140
  * McDermott & Taam (1987) McDermott, P. N., & Taam, R. E. 1987, ApJ, 318, 278
  * Numata & Lee (2010) Numata, K., & Lee, U. 2010, MNRAS, 409, 481
  * Owen et al. (1998) Owen, B. J., Lindblom, L., Cutler, C., et al. 1998, Phys. Rev. D, 58, 084020
  * Passamonti & Andersson (2011) Passamonti, A., & Andersson, N. 2011, MNRAS, 413, 47
  * Piro & Bildsten (2004) Piro, A. L., & Bildsten, L. 2004, ApJ, 603, 252
  * Piro & Bildsten (2005) —. 2005, ApJ, 619, 1054
  * Rawls et al. (2011) Rawls, M. L., Orosz, J. A., McClintock, J. E., et al. 2011, ApJ, 730, 25
  * Strobel et al. (1999) Strobel, K., Schaab, C., & Weigel, M. K. 1999, A&A, 350, 497
  * Strohmayer & Mahmoodifar (2014) Strohmayer, T., & Mahmoodifar, S. 2014, ApJ, 784, 72
  * Strohmayer & Lee (1996) Strohmayer, T. E., & Lee, U. 1996, ApJ, 467, 773
  * Strohmayer & Markwardt (2002) Strohmayer, T. E., & Markwardt, C. B. 2002, ApJ, 577, 337
  * Strohmayer et al. (1998) Strohmayer, T. E., Zhang, W., Swank, J. H., White, N. E., & Lapidus, I. 1998, ApJ, 498, L135
  * Tout et al. (1996) Tout, C. A., Pols, O. R., Eggleton, P. P., & Han, Z. 1996, MNRAS, 281, 257
  * Wood et al. (1991) Wood, K. S., Norris, J. P., Hertz, P., et al. 1991, ApJ, 379, 295
  * Yoshida & Lee (2000) Yoshida, S., & Lee, U. 2000, ApJS, 129, 353