# Chandra Observations of Six Peter Pan Disks: Diversity of X-ray-driven
Internal Photoevaporation Rates Doesn’t Explain Their Rare Longevity

Stefan Laos Homer L. Dodge Department of Physics and Astronomy, University of
Oklahoma
440 West Brooks Street, Norman, OK 73019, USA  John P. Wisniewski Homer L.
Dodge Department of Physics and Astronomy, University of Oklahoma
440 West Brooks Street, Norman, OK 73019, USA  Marc J. Kuchner NASA Goddard
Space Flight Center, Exoplanets and Stellar Astrophysics Laboratory, Code 667,
Greenbelt, MD 20771, USA  Steven M. Silverberg  MIT Kavli Institute for
Astrophysics and Space Research, 77 Massachusetts Avenue, Cambridge, MA 02139,
USA  Hans Moritz Günther MIT Kavli Institute for Astrophysics and Space
Research, 77 Massachusetts Avenue, Cambridge, MA 02139, USA  David A. Principe
MIT Kavli Institute for Astrophysics and Space Research, 77 Massachusetts
Avenue, Cambridge, MA 02139, USA  Brett Bonine  Homer L. Dodge Department of
Physics and Astronomy, University of Oklahoma
440 West Brooks Street, Norman, OK 73019, USA  Marina Kounkel Department of
Physics and Astronomy, Vanderbilt University, VU Station 1807, Nashville, TN
37235, USA  The Disk Detective Collaboration

###### Abstract

We present Chandra X-ray observations of 6 previously-identified Peter Pan
objects, rare $\sim$40 Myr systems with evidence of primordial disk retention.
We observe X-ray luminosities (0.8-3.0 keV) ranging from log $L_{x}$ $\sim$
27.7–29.1. We find that our Peter Pan sample exhibits X-ray properties similar
to that of weak-lined T-Tauri stars and do not exhibit evidence of stellar
accretion induced X-ray suppression. Our observed Peter Pan X-ray luminosities
are consistent with that measured for field dM stars of similar spectral type
and age, implying their long primordial disk lifetimes are likely not a
consequence of unusually faint X-ray host stars. Our derived X-ray
photoevaporative mass loss rates predict our systems have passed the point of
rapid gas dispersal and call into question the impact of this internal
mechanism for primordial disk dispersal around dM stars. Our qualitative
assessment of the surrounding Peter Pan environments also does not predict
unusually low levels of external photoevaporation relative to other respective
moving group members. Overall, our results suggest Peter Pan disks may be a
consequence of the low FUV flux incident on the disk in low-mass DM stars
given their relatively lower levels of accretion over the course of their pre-
main-sequence evolution.

##  1 Introduction

Large, sensitive IR surveys conducted by NASA’s Spitzer and WISE satellites
have demonstrated both that primordial protoplanetary disks of dust and gas
clear on rapid ($\lesssim$10 Myr) timescales, and that this dissipation
timescale depends on the mass of the host protostar (Carpenter et al. 2006;
Wahhaj et al. 2010; Pecaut & Mamajek 2016). Although dispersion exists between
different environments, the FEPS and c2d legacy Spitzer surveys have found
that $\sim$10 Myr is an upper limit for the lifetimes of primordial disks
around solar-type stars (Carpenter et al. 2006; Wahhaj et al. 2010; Williams &
Cieza 2011). For lower-mass K-type stars, the frequency of primordial disks is
$\sim$9% at 10 Myr and $\sim$4% by 16 Myr, indicating a disk e-folding
timescale of $\sim$4-5 Myr (Pecaut & Mamajek, 2016). These disk lifetimes not
only set the timescale for the star formation process but also necessarily
constrain the timescale for planet formation in these systems. Understanding
the detailed evolution of primordial disks around the most abundant low-mass
stars (e.g., M dwarfs) is of strong current interest, due in part to the
expected windfall of M dwarf planet science from NASA’s TESS mission. TESS has
already confirmed 18 planets having radii less than 2 R⊕ surrounding low-mass
M dwarfs (e.g. Crossfield et al. 2019, Demangeon et al. 2021).

Once depleted of gas, primordial disks evolve into second-generation debris
disks. They are composed primarily of dust, created by continuous collisional
destruction of cometary and asteroidal bodies that are remnants of the star
and planet formation process. This destruction is thought to be at least
partially stirred by the presence of planetary bodies in the system (Wyatt,
2008). These systems have depleted most of their gas and no longer exhibit
signatures of accretion. The frequency of debris can be as high as 20% for
young main sequence A-type stars, and such disks can persist for hundreds of
Myrs (Rieke et al., 2005). By contrast, Binks 2016 found zero M dwarf debris
disks in their ALLWISE study of moving groups older than 40 Myr, consistent
with the small number of candidate M dwarf debris disk systems identified by
Theissen & West 2014 (of which some were false-positive detections, Silverberg
et al. 2018). The observed low frequency of old ($\geq$40 Myr) M dwarf debris
disks may indicate that such systems are ultimately cleared by stellar wind
(Wyatt, 2008) or stellar activity (e.g. AU Mic; Grady et al. 2020) over
shorter timescales.

Over the past few years, a new population of low-mass, M-type stars, with ages
$\sim$45 Myr reliably secured via moving group membership, have been found to
exhibit substantial (LIR/L∗ $\sim$0.1) mid-IR excesses, indicative of the
presence of warm dust disks (Silverberg et al. 2016; Boucher et al. 2016;
Murphy et al. 2018; Lee et al. 2020; Silverberg et al. 2020). While such
systems were initially interpreted to be (rare) examples of the oldest dM-type
debris disk systems (Silverberg et al., 2016), their large IR excesses suggest
there should be gas in the system, and is less consistent with the
interpretation of a debris disk origin. Indeed, growing evidence shows some of
these systems exhibit clear accretion tracers, with strong (10-125 Å), broad
(200-350 km/s) Hα emission (Boucher et al. 2016; Murphy et al. 2018;
Silverberg et al. 2020) and variable Paschen $\beta$ and Brackett $\gamma$
emission (Silverberg et al., 2020) in the canonical case of WISEA
J080822.18–644357.3. With evidence of ongoing gas accretion despite their
significantly advanced ages, this set of objects have been coined as “Peter
Pan” objects, given their refusal to “grow up” (Silverberg et al., 2020).
These systems may also be the pre-cursors to notable systems like TRAPPIST-1
(Gillon et al., 2017), in which a long-surviving disk around a low-mass dM
star circularizes numerous planetary orbits into a stable configuration.

Overall, these systems currently challenge the standard paradigm for
primordial disk (and planet formation) lifetimes. Silverberg et al. 2020 argue
the simplest explanation for these systems is that they are primordial disks
that are still in the process of dissipating, implying lower than average disk
gas-mass loss rates. A key driver of primordial-disk dissipation is
photoevaporation of the disk’s gas content from high energy photons radiated
by the host star (Hollenbach et al. 1994; Clarke et al. 2001; Alexander et al.
2004; Font et al. 2004; Ercolano et al. 2009; Gorti & Hollenbach 2009; Owen et
al. 2010; Williams & Cieza 2011). Studies of the low-mass population of the TW
Hya association have shown an anti-correlation between X-ray luminosity and
disk fraction as a function of spectral type; earlier M dwarfs have higher
X-ray luminosities and a lower disk fraction than mid-late M dwarfs (Kastner
et al., 2016). The fundamental question remains: are the disks around Peter
Pan stars primordial and if so, why have they persisted around these dwarf M
stars and not others? Do Peter Pan objects represent a rare class of X-ray
faint stars?

We investigate these open questions with our new X-ray observations, which
gauge the soft (0.8-3.0 keV) X-ray photons radiating from the host star. In
this paper, we characterize the X-ray properties of identified Peter Pan
objects and estimate their expected mass loss rates from X-ray driven
photoevaporation for the first time. In Section 2, we describe the nature of
our Chandra observations as well as our utilized data reduction and spectral
extraction techniques. In Section 3, we report our derived Peter Pan X-ray
luminosities, comparing to both active, young stars and also other dwarf
M-type stars in nearby moving groups. We also compare our findings to the
X-ray luminosities used in the recent modeling of Wilhelm & Portegies Zwart
2022. In Section 4, we interpret these observations in the context of recent
modeling and proposed disk dispersal processes. Finally, in Section 5, we
conclude with a brief summary of our findings.

##  2 Observations and Data Reduction

X-ray observations (Table 1) of our Peter Pan sample were obtained with the
Advanced CCD Imaging Spectrometer (ACIS) onboard the Chandra X-ray
Observatory. All data were acquired on ACIS-S3 chip in VFAINT mode.

To extract our X-ray spectra, we analyzed the pipeline-processed data files
provided by the Chandra X-ray Center using standard science threads with CIAO
version 4.13 and calibration data from CALDB version 4.9.5 (Fruscione et al.,
2006). Spectra were extracted using circular apertures centered on the stellar
X-ray sources, with diameters of 4″ in the case of J0501 and J0808. In the
case of close binaries J0949 (a$\sim$1.5″) and J0446 (a$\sim$2.3″), we closely
inspect the binned level 2 event files by eye, adopting optimal aperture
diameters of $\sim$1.4″and $\sim$2.6″, respectively. These sizes minimized
companion contamination while also maximizing the enclosed energy fraction
given the Chandra on-axis point spread function (PSF). For each observation,
we sample the background with a circular aperture in a nearby region on the
same CCD to each source.

To derive X-ray luminosities for our Peter Pan sample, we analyze our observed
X-ray spectra (Section 2) using version 4.13 of CIAO/SHERPA (Doe et al., 2007)
and version 4.9.5 of CALDB (Fruscione et al., 2006). We note our X-ray count
data peaks around 1 keV for all objects, expected from their young, pre-main-
sequence status. For each individual observation, we group the data to 5
counts per bin and consider counts in the 0.8-3 keV range for fitting. For our
spectral model, we adopt the XSPEC optically thin plasma model vapec,
considering only a single temperature component to avoid over-fitting our low-
count data. We also add a multiplicative photoelectric absorption component
wabs to our spectral model. We find negligible intervening absorption,
expected given the close proximity of our sources ($\lesssim$100pc), and
choose to freeze its value to zero.

Similar to the analysis done in Kastner et al. (2016), we fix our plasma
abundances to the typical values determined for T Tauri stars in Taurus
(Skinner & Güdel 2013 and references therein). Relative to solar, these
abundance values are H = 1.0, He = 1.0, C = 0.45, N = 0.79, O = 0.43, Ne =
0.83, Mg = 0.26, Al = 0.50, Si = 0.31, S = 0.42, Ar = 0.55, Ca = 0.195, Fe =
0.195, and Ni = 0.195. In the case of J0446A/B and J0808, we perform
simultaneous fits over the numerous observations of each object. To gauge the
quality of our fits, we consider the Cash statistic (Cash, 1979) as
implemented in SHERPA rather than a chi-squared metric, which is not
appropriate for the Poisson-distributed low-count regime of our data. The
resulting best fit parameters derived from this analysis are given in Table 2.
From these best fits, we derive X-ray fluxes in the 0.3-8 keV energy band and
use Gaia DR3 distances (Gaia Collaboration et al., 2021) to derive our Peter
Pan X-ray luminosities ($L_{x}$), also reported in Table 2. Our best fits are
shown in Appendix Section A. Given the low count rates and correspondingly
large relative Poisson uncertainty, we note we cannot exclude variability of a
factor of two between our observations, but would have detected order-of-
magnitude flares had they occurred during our observations.

To compute overall fractional luminosities, we also compute the bolometric
luminosity of our Peter Pan sources. Following the methodology of Kastner et
al. (2016), we use known spectral types, Gaia DR3 distances, and J band data
based on the (spectral-type-dependent) J band bolometric corrections
determined by Pecaut & Mamajek (2013) assuming no reddening. For spectral
types later than M5, we linearly extrapolate from the bolometric J correction
relation. In the cases of J0446 and J0949, we assume the flux contribution
from each stellar component is equal. These values are reported in Table 2.

Table 1: Peter Pan X-ray Observations Object | Obs. Date (UTC) | Obs. ID | Exposure (ks)
---|---|---|---
WISEA J094900.65–713803.1 A/B | 2019-12-13 15:34:31 | 22305 | 17
2MASS J05010082–4337102 | 2019-12-22 5:15:02 | 22306 | 19
WISEA J044634.16–262756.1 A/B | 2020-05-20 4:48:40 | 22304 | 22
| 2020-05-21 8:57:52 | 23255 | 10
| 2020-06-02 4:27:50 | 23256 | 23
WISEA J080822.18–644357.3 | 2021-07-06 14:29:42 | 22303 | 36
| 2021-07-07 6:24:58 | 24756 | 35
Table 2: X-ray Spectral Analysis: Best-fit Parameters Source | Fx | kT | norm | Reduced CstataaThe Cstat procedure is the XSPEC111http://heasarc.gsfc.nasa.gov/xanadu/xspec. implementation of the Cash statistic (Cash, 1979), which also incorporates fitting a model to the background. The reduced cstat values correspond to the observed statistic divided by the degrees of freedom, operating as an approximate goodness of fit metric with good fits to the data approaching a nominal value of 1. | Lx | Lbol | Spectral TypebbSpectral types are derived from the analysis of Silverberg et al. 2020 and references therein. Ages are estimated from the isochronal fitting of Bell et al. 2015. | AgebbSpectral types are derived from the analysis of Silverberg et al. 2020 and references therein. Ages are estimated from the isochronal fitting of Bell et al. 2015.
---|---|---|---|---|---|---|---|---
| erg/s/cm2 | keV |  |  | erg/s | erg/s |  | Myr
J0949A | 1.48$\pm$0.4E-13 | 0.28$\pm$0.02 | 3.0E-4$\pm$7.5E-5 | 1.69 | 1.1E+29 | 2.18E+32 | M4 | 45
J0949B | 6.55$\pm$0.9E-14 | 0.73$\pm$0.13 | 8.4E-5$\pm$1.5E-5 | 1.39 | 4.8E+28 | 1.92E+32 | M5 | 45
J0501 | 2.14$\pm$0.6E-14 | 0.63$\pm$0.28 | 2.0E-5$\pm$9.2E-6 | 1.49 | 5.9E+27 | 2.60E+31 | M4.5 | 42
J0446A | 4.39$\pm$0.4E-14 | 0.95$\pm$0.09 | 3.8E-5$\pm$3.4E-6 | 1.05 | 3.6E+28 | 5.24E+31 | M6 | 42
J0446B | 3.03$\pm$0.9E-14 | 0.50$\pm$0.10 | 2.7E-5$\pm$6.3E-6 | 1.65 | 2.5E+28 | 5.19E+31 | M6 | 42
J0808 | 1.21$\pm$0.4E-14 | 0.46$\pm$0.12 | 1.1E-5$\pm$3.7E-6 | 0.77 | 1.5E+28 | 3.08E+31 | M5 | 45

##  3 Results

With our derived luminosities, we explore the X-ray properties of our 6 Peter
Pan targets for the first time. In Section 3.1, we consider our X-ray
luminosities in the context of stellar accretion. In Section 3.2, we compare
our derived values to sensitive X-ray observations of our M star analog
control sample (Appendix Section B). We also compare our luminosities to that
utilized in the prescriptions of recent modeling in Section 3.3.

###  3.1 Accretion vs. Lx

We begin by comparing our sample with the large sample of classical and weak-
lined T-Tauri stars (CTTS, WTTS) in Taurus from Telleschi et al. 2007. The
sensitive XMM-Newton observations of this young (average age of $\sim$2.4
Myr), active sample allow us to further gauge the impact of active stellar
accretion on $L_{x}$.

Probing even lower bolometric luminosities, our observations find a large
diversity in $L_{x}$ across the Peter Pan sample (Figure 1), spanning a range
of $\sim$1.5 dex. Their overall appreciable X-ray luminosities do not show
evidence for the significant ($\sim$2x) X-ray suppression seen in the
(actively accreting) CTTS subsample (Telleschi et al., 2007). We instead find
the $L_{x}$/$L_{bol}$ values of the Peter Pan sources more closely follow the
trend observed in the WTTS, consistent with the similarly weak level of Peter
Pan accretion ($\sim$10-11-10-9 M$\odot$ yr-1) derived from recent H-alpha
observations Silverberg et al. 2020.

We also compare the H-alpha emission equivalent widths from the measurements
of Silverberg et al. 2020 with our observed X-ray luminosities. We do not find
evidence of a strong correlation (computing a Pearson correlation coefficient
value of -0.45) within the limitations of our small sample. Interestingly, we
note the strongest accretor (J0949A) is also our brightest X-ray object,
counter to the expectation from stellar accretion induced X-ray suppression.
Collectively, these results do not find evidence of significant X-ray
suppression in the Peter Pan sample as seen in the active CTTS, likely due to
their relatively low level of accretion.

Figure 1: Comparison of observed X-ray luminosities between our Peter Pan
sample (starred points) and the T Tauri stars from Telleschi et al. 2007.
Upper limits are indicated with arrows. Diagonal black lines represent
constant X-ray to bolometric luminosity fractions (Lx/Lbol) of 10-3, 10-4,
10-5. Linear regressions for the WTTS (green) and CTTS (purple) subsamples are
reproduced from Telleschi et al. 2007, shown as dotted lines. We find the
Peter Pan sources more closely resemble the Lx/Lbol behavior exhibited by the
WTTS (green) relative to the X-ray fainter CTTS (purple).

###  3.2 M Star Control Sample Comparison

Our comparison with the T Tauri stars from (Telleschi et al., 2007) in Section
3.1 mainly considers relatively young stars ($\sim$1-10 Myr). Given the
suspected advanced ages of the Peter Pan sources ($\sim$45 Myr, Silverberg et
al. 2020), we construct a new comparison sample, closely matching suspected
spectral types and ages with similarly sensitive X-ray observations. We
describe our control sample criteria in Appendix Section B.

In Figure 2, we compare both X-ray (left panel) and fractional X-ray (right
panel) luminosities between our Peter Pan sources and this control sample.
Overall, we find the values observed are consistent with one another. We
therefore do not find evidence that these known Peter Pan sources are over- or
under-luminous in X-rays relative to similarly evolved M stars. This finding
may argue against the previously thought rarity of disk longevity out to Peter
Pan ages for mid-M types. We explore this possibility further in Section 4.

Figure 2: Observed X-ray (left) and fractional X-ray (right) luminosities as a
function of spectral type between the Peter Pan sample and M star control
sample (i.e., similar ages and spectral types of our Peter Pan sample).
Circular points represent spectroscopically verified spectral types while
crosses represent spectral types estimated via colors. Overall, we find our
derived Peter Pan luminosities to be consistent with that observed for other
moving group M stars.

###  3.3 Predicted vs. Observed X-ray Luminosity Evolution in dM Stars

Our X-ray observations offer a set of measurements to benchmark the
characteristic X-ray luminosity for Peter Pan disk systems and assess
previously used predictions. Recent modeling has investigated the survival of
primordial circumstellar disks out to the advanced ages of 45 Myr (Wilhelm &
Portegies Zwart 2022), assuming disk winds from strong stellar X-rays
constitute the primary internal source of disk dispersal. Detailed simulations
of this mass loss procedure (Owen et al. 2012; Picogna et al. 2019) have found
a power-law dependence on the X-ray luminosity of the host star (which in turn
is dependent on stellar mass).

To estimate the expected Peter Pan X-ray luminosity, Wilhelm & Portegies Zwart
2022 consider the characteristic (quiescent) X-ray luminosity-mass relation
derived by Flaccomio et al. 2012 for young stars observed by the Chandra Orion
Ultra Deep Project (COUP, Getman et al. 2005), which is nearly complete
($\sim$95%, Preibisch et al. 2005) down to the stellar mass limit. These X-ray
luminosities steadily decrease throughout pre-main-sequence evolution,
estimated to follow a time dependence of L${}_{x}\propto t^{-2/5}$ (Gregory et
al. 2016; Johnstone et al. 2021). Together, these results give the functional
form given in Equation 1, where a and b are constants 30.0 and 1.87
respectively (details noted in Section 2.2 of Wilhelm & Portegies Zwart 2022).

$L_{x}(t)=10^{a+(b*\text{log}[M_{*}/M\textsubscript{$\odot$}])}\left(\frac{t}{1\text{Myr}}\right)^{-2/5}$
(1)

In Figure 3, we compare our observed Peter Pan X-ray luminosities with the
predictions from Equation 1, assuming main sequence masses and the nominal
ages estimated for the Peter Pan systems in Silverberg et al. 2020. Overall,
we find this prescription roughly estimates the locus of our measurements but
in some cases differs by up to $\sim$0.5 dex in the cases of our brightest
(J0949A) and faintest (J0501) Peter Pan source. We note however that a more
robust comparison is difficult, given the large uncertainties on our derived
x-ray luminosities and the small size of our sample. This is further
complicated by the large intrinsic scatter in observed stellar x-ray
luminosities, varying by factor of $\sim$2-3 as a consequence of variations in
surface magnetic activity, magnetic cycles, or flares (Preibisch et al. 2005,
Johnstone et al. 2021). Although we find no evidence of flares in our x-ray
light curve data, it remains unclear the extent to which Peter Pan x-ray
luminosities are being over- or under- predicted in recent modeling. We
discuss the implications of this result and investigate their expected mass
loss rates further in Section 4.1.

Figure 3: Comparison of observed Lx against the predicted Lx values used in
the recent modeling conducted by Wilhelm & Portegies Zwart 2022 (Section 3.3).
The dotted line shown represents the 1-1 relation. 1$\sigma$ error bars are
shown in light red.

##  4 Discussion

Primordial disks are subject to mass loss through a variety of mechanisms,
both internal and external to the system (Armitage 2011; Williams & Cieza
2011). In the inner disk, these mechanisms include both the accretion of disk
material onto the stellar surface (Hartmann et al., 2016) as well as disk
outflows driven by the high energy flux from the central star, a process known
as internal photoevaporation (Hollenbach et al. 1994; Clarke et al. 2001;
Alexander et al. 2004; Font et al. 2004; Ercolano et al. 2009; Gorti &
Hollenbach 2009; Owen et al. 2010; Williams & Cieza 2011). Photoevaporation is
also capable of externally triggering mass loss in the outer disk as a
consequence of the FUV radiation from nearby massive stars (Störzer &
Hollenbach 1999; Scally & Clarke 2001; Adams et al. 2004; Concha-Ramírez et
al. 2019; Winter et al. 2019).

The relative importance of these internal and external channels driving disk
dissipation are suspected to depend on stellar mass. Higher mass
($\gtrsim$1M$\odot$) stars have similarly higher levels of accretion
(Alexander & Armitage 2006; Alcalá et al. 2014) and suspected internally
driven photoevaporation (Owen et al. 2012). These mechanisms likely dominate
in this case, given primordial disks around solar-type stars dissipate by
$\sim$10 Myr regardless of their diverse surrounding environments (Carpenter
et al. 2006; Wahhaj et al. 2010; Pecaut & Mamajek 2016). In the lower mass
($\lesssim$0.3M$\odot$) regime of the Peter Pan stars, the impact of external
photoevaporation becomes more comparable even in low external radiation fields
(Haworth et al., 2018) given their lower internal accretion. Overall, the
relative impact of internal and external mechanisms remain debated for the
lowest-mass stars, with these respective class of mechanisms likely dominating
disk dispersal at different points of pre-main-sequence evolution (Wilhelm &
Portegies Zwart 2022). We discuss the implications of our X-ray results in
regards to the described internal and external mechanisms in Section 4.1 and
4.2, respectively.

###  4.1 The Impact of Internal Disk Dispersal Mechanisms on Peter Pan Disk
Longevity

The overall effect of internal photoevaporation on disk dispersal has remained
largely unconstrained, as investigations continue to debate the overall,
relative importance of each stellar flux energy regime (e.g. far-ultraviolet,
extreme-ultraviolet, and X-ray) on resultant mass-loss rates and disk
lifetimes. Recent work has argued the impact of EUV radiation is limited,
given its absorption over small column densities that does not allow for
penetration into the high-density disk midplane (Owen et al. 2010; Picogna et
al. 2019). In general, the wind mass-loss predictions from EUV-dominated
photoevaporation studies (Hollenbach et al. 1994; Tanaka et al. 2013) have
also been lower by an order of magnitude or more than that from FUV and X-ray
investigations (Ercolano et al. 2009; Gorti & Hollenbach 2009; Owen et al.
2010, 2012; Gorti et al. 2015). Given this and the nature of our observations,
we choose to mainly interpret our results in the context of current
predictions from X-ray driven disk photoevaporation modeling.

Wind mass-loss rates from X-ray driven photoevaporation have been studied in
growing detail, with a library of models covering the observed parameter space
of both stellar X-ray properties as well as disk metallicity (Picogna et al.
2019; Wölfer et al. 2019). Ercolano et al. 2009 found that the stellar
luminosity in the soft (0.1-1 keV) X-ray band drives the bulk of the
photoevaporative wind, as harder X-rays are unable to provide enough disk
heating to unbind their gaseous components despite deeper penetration (Wang &
Goodman 2017; Nakatani et al. 2018). Ercolano et al. 2021 derive a functional
form (Equation 2) for the total expected mass loss from X-ray photoevaporation
(given a soft X-ray luminosity down to log L${}_{x}\sim$28.5) and find their
estimates are in agreement with previous modeling (Owen et al. 2012; Picogna
et al. 2019). Prefactor constants have values of aS = -1.947 x 1017, bS =
-1.572 x 10-4, cS = -2.866 x 10-1, and dS = -6.694, respectively.

$\text{log}\dot{M}_{W}(L_{X,\text{soft}})=a_{S}~{}\text{exp}\left(\frac{(\text{ln}(\text{log}L_{X,\text{soft}})-b_{S})^{2}}{c_{S}}\right)+d_{S}$
(2)

Figure 4: Comparison of the expected mass loss rates between accretion and
X-ray driven photoevaporative disk winds for our Peter Pan sample. Values for
the former are reported in Silverberg et al. 2020 while values for the latter
are computed from the functional form (Equation 2) derived in Ercolano et al.
2021. We indicate the 1-1 relation with a black dashed line.

Using our best fit models in the 0.8-3.0 keV band, we extrapolate to determine
the 0.1-1.0 keV X-ray luminosities of our Peter Pan sample and calculate their
wind mass-loss rates from Equation 2 (as detailed in Ercolano et al. 2021). We
report these values in Table 3, conservatively estimating an uncertainty of
$\sim$1 dex. We find our observed X-ray luminosities predict appreciable wind
mass-loss rates that also span a large range (over 2 dex) across our entire
sample. We note these values are effectively lower limits given their X-ray
luminosities have decreased throughout their pre-main-sequence evolution.

These derived wind loss values track the flow of inner disk material ejected
into the outer regions of the disk, in contrast to the on-going stellar
accretion, which funnels inner disk material to the stellar photosphere. As
the mass accretion rate decreases and becomes comparable to the wind mass loss
rate, the opening of a gap in the inner disk is expected, clearing the inner
disk of material on rapid timescales ($\lesssim$1 Myr, Williams & Cieza 2011).
We compare our X-ray photoevaporative mass loss rates to the Peter Pan
accretion rates (derived in Silverberg et al. 2020 from H-alpha measurements)
to gauge this critical moment of inner disk evolution in Figure 4. Overall, we
find X-ray photoevaporative mass loss rates that are larger than those derived
from accretion. Their relative difference could be a consequence of numerous
interpretations, which we discuss below.

The relatively low mass loss rates derived from accretion are not surprising
given the Peter Pans represent a sample of fairly evolved low-mass stars.
Wilhelm & Portegies Zwart 2022 have shown that lower disk viscosities, which
slow the overall accretion, can allow for longer primordial disk lifetimes
(upwards of $\sim$50 Myr for estimated disk properties around stars of
$\lesssim$0.5M$\odot$). However, assuming the rates derived in Equation 2 are
correct, these results suggest gaps are expected for all of our Peter Pan
systems. The two-component SED fits found for J0808 (Murphy et al. 2018;
Flaherty et al. 2019) could be evidence of a gap in a Peter Pan disk
supporting this interpretation.

The stronger outliers observed in J0446A/B and J09494A/B are more difficult to
explain. We note the greater relative difference in the mass loss rates for
J0446A and J0446B could be evidence of H-alpha emission sourced from stellar
activity and not accretion. These systems exhibited H-alpha equivalent widths
(10-17Å) and velocity widths (210-240 km/s) that bordered the criteria for an
accretion interpretation (Jayawardhana et al. 2003; Fang et al. 2009). While
our X-ray results may argue against the Peter Pan status of these systems,
this explanation is unable to resolve the tension with J0949A, which is much
more clearly accreting given its strong (110Å) and broad ($\sim$370 km/s)
emission (Silverberg et al., 2020).

X-ray photoevaporative models predict this disk gap phase occurs at an age of
$\sim$10-20 Myr (Ercolano et al. 2021; Picogna et al. 2021) and constitutes
the last $\sim$25% of the disk lifetime (for a stellar X-ray luminosity of
1029 erg s-1), leading to a quick dispersal of the remaining disk gas in only
a few Myr (Ercolano et al. 2021). These short remaining primordial disk
lifetimes would be in tension overall with the suspected Peter Pan ages
($\sim$40 Myr). Given it is also unlikely to have caught these systems at a
point at which X-ray photoevaporative rates have spiked, we argue the
$\dot{M}_{W}$–Lx relation found by Ercolano et al. 2021 is likely inaccurate
for our Peter Pan sample. The discrepancy may be explained by the fact that
our observed Peter Pan luminosity range (log $L_{x}$ $\sim$ 27.7–29.1) extends
much lower than the range tested in the models of Ercolano et al. 2021 (log
$L_{x}$ $\sim$ 29–31).

Alternatively, our derived X-ray photoevaporative wind mass loss rates could
be overestimated in the case of mechanisms that shield the inner disk from the
stellar X-ray radiation. It has been hypothesized that hard X-ray photons
could be screened from interacting with the disk via a molecular magneto-
hydrodynamic inner disk wind, driven by the accretion-powered active stellar
magnetic field (Pascucci et al., 2020). This mechanism, however, is mainly
predicted for strong classical accretors (10-8 Msun yr-1) and is likely
inconsistent with early Peter Pan evolution, given their relatively low
masses.

Lastly, we note recent findings have suggested younger ages for the Carina and
Columba moving groups, ranging from 13 (Booth et al., 2021) to 22 (Schneider
et al., 2019) Myr. These ages could reduce the tension in their current
classification, given a disk e-folding timescale of $\sim$4-5 Myr (as observed
for K-type stars, e.g. Pecaut & Mamajek 2016) could explain a small fraction
of $\sim$20 Myr primordial disks. However, these results are in disagreement
with the lithium depletion boundary analysis of Murphy et al. 2018, which
indicate a suspected age of $\sim$40 Myr for J0808. These younger moving group
ages would also not be able to explain the Peter Pan disk candidates
identified in other moving groups, such as 2MASS J02265658–5327032 in the
$\sim$45 Myr Tuc-Hor (Silverberg et al., 2020) moving group and 2MASS
J15460752-6258042 in the $\sim$55 Myr Argus moving group (Lee et al., 2020).
Continued efforts to determine more exact ages for these moving groups are
needed.

Overall, our findings question the impact of X-ray photoevaporative disk winds
as the primary source of disk dispersal in dM stars. We note photoevaporative
winds driven by the host stellar FUV radiation cannot be ruled out from our
observations. It has been suspected a large fraction of FUV photons are
sourced from the accretion shocks, with a broad range of FUV luminosities
($10^{-6}$L$\odot$ $\lesssim$ LFUV $\lesssim$ L$\odot$) observed in active
classical TTS (Gullbring et al. 1998; Yang et al. 2012). Therefore, Peter Pan
disks may be a consequence of the low FUV flux incident on the disk of low-
mass dM stars given their relatively low levels of accretion over the course
of their pre-main-sequence evolution. Future FUV observations of low-mass dM
stars will be critical to confirm this behavior. If FUV photoevaporation is
indeed critical for dM stars, however, it remains unclear why more Peter Pan
disks have yet to be discovered.

Table 3: Mass Loss Rate Estimations Source | Lx (0.1-1 keV) | log (Macc)aaEstimated in Silverberg et al. 2020 from H-alpha observations. | log (Mw)bbComputed using the $\dot{M}_{W}$–Lx relation derived by Ercolano et al. 2021 (Equation 2).
---|---|---|---
| erg/s | Msun / yr | Msun / yr
J0949A | 9.77E+28 | -9.3 | -7.8
J0949B | 2.81E+28 | -9.9 | -8.2
J0501 | 4.15E+27 | -10.8 | -10.2
J0446A | 2.12E+28 | -10.9 | -8.6
J0446B | 1.88E+28 | -10.6 | -8.8
J0808 | 1.17E+28 | -9.75 | -9.2

###  4.2 The Impact of External Photoevaporation Mechanisms on Peter Pan Disk
Longevity

The intense FUV and EUV fields radiating from the rare, massive O and B stars
in a dense star forming region serve as another source of disk dispersal by
driving a photoevaporative wind in the outer circumstellar disk. The overall
impact and efficiency of this external mechanism on primordial disk longevity
is still currently debated. Previous observational studies of young (1-5 Myr)
star forming regions have found evidence of a correlation between the
projected distance to the most luminous ionizing stellar source and the
suspected mass of nearby disks in some cases (e.g. Orion Nebula Cluster: Mann
et al. 2014, Eisner et al. 2018; $\sigma$ Orionis: Ansdell et al. 2017, van
Terwisga et al. 2019) but not others (e.g. NGC 2024: Mann et al. 2015). Recent
N-body simulations have not found any evidence of a correlation (Nicholson et
al. 2019) even after accounting for projection effects (Parker et al. 2021).
Similarly, disk lifetimes were not found to be a function of projected
distance to the cluster core in IC 1396 (Silverberg et al., 2021). Even so,
these investigations do find that the extreme case of disks in close
($\sim$0.5 pc) proximity to massive stars are effectively dissipated, with few
disks surviving longer than a few Myr on average.

The evidence of $\sim$45 Myr old primordial disk systems in our Peter Pan
sample inherently implies significantly low levels of photoevaporation from
external radiation for prolonged durations. The investigations of Coleman &
Haworth 2020 find primordial disks around low mass stars could survive up to
$\sim$50 Myr in the presence of low external radiation fields ($<$10 Habing
fields (1.6 x 10-2 erg cm-2 s-1), corresponding to an external
photoevaporative mass loss rate $\leq$10-9 M$\odot$/yr). A direct quantitative
estimation of the local radiation experienced by our Peter Pan systems is not
feasible. We note qualitative assessments of the external radiation field,
based off of the current 3d positions of the Peter Pan sample relative to
other moving group members, are also difficult to interpret. The spatial
stellar distribution is not informative of the respective moving group
dynamics soon after formation, times at which the stellar density was higher
and close encounters more probable.

Although we are unable to probe these early times ($\lesssim$2 Myr), we note
the spatial and kinematic properties of Carina, Columba, and Tuc-Hor are
somewhat suggestive of unique dynamical evolution, given their fairly large
spatial extents ($\sim$10-15 pc) and particularly low velocity dispersions
($\lesssim$1 km/s) relative to other known moving groups (Gagné et al., 2018).
Recent evidence has also suggested that nearby moving groups of stars may be
spatially much more extended than previously thought. The high-precision
stellar kinematics offered by Gaia DR2 have revealed numerous co-moving
stellar streams throughout the galaxy at distances ranging 80-1000 pc (Kounkel
& Covey 2019). Some of these “Theia” strings have been identified as
relatively young ($\lesssim$100 Myr) and are suspected to be primordial,
reflecting the shape of the molecular cloud from which they have formed. In
some cases, these streams appear to be interconnected with nearby moving
groups, exhibiting similar estimated isochoronal ages and Gaia DR2 UVW space
velocities (Gagné et al. 2021). We test for Theia string membership among our
Peter Pan sample and note a consistent match with Theia 113 for J0808.

Although the overall dynamical evolution of moving groups and respective Theia
strings is still not well understood, $\sim$2-3 crossing times are still
predicted for the stars in these moving groups. The expectation is then that
all moving group members are equally likely to have experienced nearby
proximity to massive members, in contrast with an interpretation of strict
isolation resulting in low external photoevaporation. The lack of confirmed
massive members as assessed by the BANYAN algorithm (Gagné & Faherty 2018) in
the moving groups of Carina, Columba, and Tuc-Hor (1 B9 star and no O stars)
also argue against the overall impact of external photoevaporation.

If external photoevaporation is most important for disk dispersal in dM stars,
we would expect the dM stars nearest the Peter Pan sources to experience
similar levels of incident external radiation, and thus also be candidates for
harboring long-lived primordial disks. We refer to our list of mid-M stars
with high-probability moving group membership (J. Gagné, private
communication) to explore this sub-sample. We find the majority of these
candidates did not pass the original criteria of the Disk Detective search
(Kuchner et al., 2016) as their disk harboring status could not be determined
due to low-SNR WISE observations. We cross-reference this sub-sample with Gaia
DR 3 and report the 5 nearest mid-M stars to the canonical Peter Pan target,
J0808, in Table 4. Follow-up near-IR observations of these low-mass stars
would help test the extent to which the long-lived disks in our Peter Pan
sample could be explained by unusual local, external radiation fields.

Table 4: Nearest Carina Mid-Ms to J0808 Source | RA | DEC | SpT | 3D SeparationaaEstimated from Gaia DR 3 positions and parallaxes.
---|---|---|---|---
| deg | deg |  | pc
2MASS J08111195-6656400 | 122.79962328091591 | -66.94434486171754 | M4.0 | 8.3
2MASS J08441995-6158424 | 131.08297932821128 | -61.97836028331041 | M4.9 | 13.2
2MASS J08122535-6852061 | 123.105464898878 | -68.8682618587011 | M7.0 | 15.5
2MASS J07151705-6555486 | 108.82103844623693 | -65.93003833077316 | M4.3 | 17.0
2MASS J07550342-6717478 | 118.76411975737 | -67.2964623058796 | M4.7 | 19.5

Overall, we do not find evidence suggestive of lower external photoevaporation
experienced by our Peter Pan sample relative to other moving group members.
However, we acknowledge that many of the aspects of both early moving group
dynamic and the efficiency of external photoevaporation in the regime of low
stellar mass, remain unconstrained. If Peter Pan disk longevity is a
consequence of an usually low external FUV radiation field, we predict Peter
Pan systems would likely require unusually high amounts of self-shielding
during their early evolution to survive to their current suspected ages of 45
Myr and explain their overall rarity.

##  5 Summary and Conclusion

We present new X-ray observations for our sample of 6 recently identified
Peter Pan stars, systems with strong evidence of harboring primordial disks at
the advanced ages of $\sim$45 Myr. Our results in Section 3 reveal the X-ray
characteristics of these sources for the first time. We summarize our main
results below.

  * •

We observe Peter Pan X-ray emission similar to that observed in young, weak-
lined TTS (Figure 1). Their derived X-ray luminosities have an overall large
dispersion, spanning roughly 1.5 dex (log Lx = 27.7-29.1). We do not find
evidence of their X-ray luminosities correlating with their H-alpha equivalent
widths, arguing against accretion-driven X-ray suppression being responsible
for their extended disk lifetimes.

  * •

Our derived Peter Pan X-ray luminosities are consistent with that measured for
field dM stars of similar spectral type and age (Figure 2). We conclude Peter
Pan disk lifetimes are likely not a consequence of central stars with lower
Lx/Lbol.

  * •

In some cases, our observed X-ray luminosities differ from those used in the
recent modeling of Wilhelm & Portegies Zwart 2022 by up to 0.5 dex (Figure 3).
It remains unclear the extent to which Peter Pan x-ray luminosities are being
over or under predicted, given our large measurement uncertainties, overall
small sample size, and the large intrinsic scatter in observed x-ray
luminosities.

  * •

Our derived X-ray photoevaporative mass loss rates (Section 4.1) are in most
cases much larger than that estimated by previous H-alpha measurements. These
measurements may suggest our Peter Pan systems are in the final phase of their
overall disk lifetime, predicted to occur at $\sim$10-20 Myr (Ercolano et al.
2021, Picogna et al. 2021). Given the disparity between this timescale and the
suspected ages of Peter Pan systems ($\sim$45 Myr), we argue the wind mass
loss relation of Ercolano et al. 2021 is inaccurate for our Peter Pan sample.
Overall, these findings call into question the impact of X-ray
photoevaporation on disk dispersal for dM stars and may suggest the importance
of FUV-driven disk dispersal mechanisms. Given the size of our sample,
additional FUV and X-ray observations of low-mass dM stars are needed to help
confirm this result.

  * •

Our qualitative assessment of the surrounding Peter Pan environments (Section
4.2) does not predict unusually low levels of external photoevaporation
relative to other respective moving group members. However, given this
assessment does not inform of early moving group dynamics, the overall impact
of external photoevaporation remains unclear. If Peter Pan disk longevity is a
consequence of an usually low external FUV radiation field, we predict that
Peter Pan systems would likely have to have experienced unusual amounts of
self-shielding to explain their current rarity.

This research has made use of the SIMBAD database and the VizieR catalogue
access tool, operated at CDS, Strasbourg, France. This research has also made
use of NASA’s Astrophysics Data System and ds9, a tool for data visualization
supported by the Chandra X-ray Science Center (CXC) and the High Energy
Astrophysics Science Archive Center (HEASARC) with support from the JWST
Mission office at the Space Telescope Science Institute for 3D visualization.
Multiple Python libraries aided the analysis of our data including matplotlib,
a Python library for publication quality graphics (Hunter, 2007), SciPy
(Virtanen et al., 2020), NumPy (Harris et al., 2020), and Astropy, a
community-developed core Python package for Astronomy (Astropy Collaboration
et al., 2013, 2018). IRAF is distributed by the National Optical Astronomy
Observatory, which is operated by the Association of Universities for Research
in Astronomy (AURA) under cooperative agreement with the National Science
Foundation (Tody, 1993). These acknowledgements were compiled using the
Astronomy Acknowledgement Generator. The scientific results reported in this
article are based on observations made by the Chandra X-ray Observatory. This
research has made use of software provided by the Chandra X-ray Center (CXC)
in the application packages CIAO (Fruscione et al., 2006) and Sherpa (Doe et
al., 2007). CXO (ACIS)

## References

  * Adams et al. (2004) Adams, F. C., Hollenbach, D., Laughlin, G., & Gorti, U. 2004, ApJ, 611, 360, doi: 10.1086/421989
  * Alcalá et al. (2014) Alcalá, J. M., Natta, A., Manara, C. F., et al. 2014, A&A, 561, A2, doi: 10.1051/0004-6361/201322254
  * Alexander & Armitage (2006) Alexander, R. D., & Armitage, P. J. 2006, ApJ, 639, L83, doi: 10.1086/503030
  * Alexander et al. (2004) Alexander, R. D., Clarke, C. J., & Pringle, J. E. 2004, MNRAS, 354, 71, doi: 10.1111/j.1365-2966.2004.08161.x
  * Ansdell et al. (2017) Ansdell, M., Williams, J. P., Manara, C. F., et al. 2017, AJ, 153, 240, doi: 10.3847/1538-3881/aa69c0
  * Armitage (2011) Armitage, P. J. 2011, ARA&A, 49, 195, doi: 10.1146/annurev-astro-081710-102521
  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33, doi: 10.1051/0004-6361/201322068
  * Astropy Collaboration et al. (2018) Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123, doi: 10.3847/1538-3881/aabc4f
  * Bell et al. (2015) Bell, C. P. M., Mamajek, E. E., & Naylor, T. 2015, MNRAS, 454, 593, doi: 10.1093/mnras/stv1981
  * Binks (2016) Binks, A. 2016, in Young Stars & Planets Near the Sun, ed. J. H. Kastner, B. Stelzer, & S. A. Metchev, Vol. 314, 159–162, doi: 10.1017/S1743921315006158
  * Booth et al. (2021) Booth, M., del Burgo, C., & Hambaryan, V. V. 2021, MNRAS, 500, 5552, doi: 10.1093/mnras/staa3631
  * Boucher et al. (2016) Boucher, A., Lafrenière, D., Gagné, J., et al. 2016, ApJ, 832, 50, doi: 10.3847/0004-637X/832/1/50
  * Carpenter et al. (2006) Carpenter, J. M., Mamajek, E. E., Hillenbrand, L. A., & Meyer, M. R. 2006, ApJ, 651, L49, doi: 10.1086/509121
  * Cash (1979) Cash, W. 1979, ApJ, 228, 939, doi: 10.1086/156922
  * Clarke et al. (2001) Clarke, C. J., Gendrin, A., & Sotomayor, M. 2001, MNRAS, 328, 485, doi: 10.1046/j.1365-8711.2001.04891.x
  * Coleman & Haworth (2020) Coleman, G. A. L., & Haworth, T. J. 2020, MNRAS, 496, L111, doi: 10.1093/mnrasl/slaa098
  * Concha-Ramírez et al. (2019) Concha-Ramírez, F., Wilhelm, M. J. C., Portegies Zwart, S., & Haworth, T. J. 2019, MNRAS, 490, 5678, doi: 10.1093/mnras/stz2973
  * Crossfield et al. (2019) Crossfield, I. J. M., Waalkes, W., Newton, E. R., et al. 2019, ApJ, 883, L16, doi: 10.3847/2041-8213/ab3d30
  * Demangeon et al. (2021) Demangeon, O. D. S., Zapatero Osorio, M. R., Alibert, Y., et al. 2021, A&A, 653, A41, doi: 10.1051/0004-6361/202140728
  * Doe et al. (2007) Doe, S., Nguyen, D., Stawarz, C., et al. 2007, in Astronomical Society of the Pacific Conference Series, Vol. 376, Astronomical Data Analysis Software and Systems XVI, ed. R. A. Shaw, F. Hill, & D. J. Bell, 543
  * Eisner et al. (2018) Eisner, J. A., Arce, H. G., Ballering, N. P., et al. 2018, ApJ, 860, 77, doi: 10.3847/1538-4357/aac3e2
  * Ercolano et al. (2009) Ercolano, B., Clarke, C. J., & Drake, J. J. 2009, ApJ, 699, 1639, doi: 10.1088/0004-637X/699/2/1639
  * Ercolano et al. (2021) Ercolano, B., Picogna, G., Monsch, K., Drake, J. J., & Preibisch, T. 2021, MNRAS, 508, 1675, doi: 10.1093/mnras/stab2590
  * Fang et al. (2009) Fang, M., van Boekel, R., Wang, W., et al. 2009, A&A, 504, 461, doi: 10.1051/0004-6361/200912468
  * Flaccomio et al. (2012) Flaccomio, E., Micela, G., & Sciortino, S. 2012, A&A, 548, A85, doi: 10.1051/0004-6361/201219362
  * Flaherty et al. (2019) Flaherty, K., Hughes, A. M., Mamajek, E. E., & Murphy, S. J. 2019, ApJ, 872, 92, doi: 10.3847/1538-4357/aaf794
  * Font et al. (2004) Font, A. S., McCarthy, I. G., Johnstone, D., & Ballantyne, D. R. 2004, ApJ, 607, 890, doi: 10.1086/383518
  * Fruscione et al. (2006) Fruscione, A., McDowell, J. C., Allen, G. E., et al. 2006, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 6270, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, ed. D. R. Silva & R. E. Doxsey, 62701V, doi: 10.1117/12.671760
  * Gagné & Faherty (2018) Gagné, J., & Faherty, J. K. 2018, ApJ, 862, 138, doi: 10.3847/1538-4357/aaca2e
  * Gagné et al. (2021) Gagné, J., Faherty, J. K., Moranta, L., & Popinchalk, M. 2021, ApJ, 915, L29, doi: 10.3847/2041-8213/ac0e9a
  * Gagné et al. (2018) Gagné, J., Mamajek, E. E., Malo, L., et al. 2018, ApJ, 856, 23, doi: 10.3847/1538-4357/aaae09
  * Gaia Collaboration et al. (2021) Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2021, A&A, 649, A1, doi: 10.1051/0004-6361/202039657
  * Getman et al. (2005) Getman, K. V., Flaccomio, E., Broos, P. S., et al. 2005, ApJS, 160, 319, doi: 10.1086/432092
  * Gillon et al. (2017) Gillon, M., Triaud, A. H. M. J., Demory, B.-O., et al. 2017, Nature, 542, 456, doi: 10.1038/nature21360
  * Gorti & Hollenbach (2009) Gorti, U., & Hollenbach, D. 2009, ApJ, 690, 1539, doi: 10.1088/0004-637X/690/2/1539
  * Gorti et al. (2015) Gorti, U., Hollenbach, D., & Dullemond, C. P. 2015, ApJ, 804, 29, doi: 10.1088/0004-637X/804/1/29
  * Grady et al. (2020) Grady, C. A., Wisniewski, J. P., Schneider, G., et al. 2020, ApJ, 889, L21, doi: 10.3847/2041-8213/ab65bb
  * Gregory et al. (2016) Gregory, S. G., Adams, F. C., & Davies, C. L. 2016, MNRAS, 457, 3836, doi: 10.1093/mnras/stw259
  * Gullbring et al. (1998) Gullbring, E., Hartmann, L., Briceño, C., & Calvet, N. 1998, ApJ, 492, 323, doi: 10.1086/305032
  * Harris et al. (2020) Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357, doi: 10.1038/s41586-020-2649-2
  * Hartmann et al. (2016) Hartmann, L., Herczeg, G., & Calvet, N. 2016, ARA&A, 54, 135, doi: 10.1146/annurev-astro-081915-023347
  * Haworth et al. (2018) Haworth, T. J., Clarke, C. J., Rahman, W., Winter, A. J., & Facchini, S. 2018, MNRAS, 481, 452, doi: 10.1093/mnras/sty2323
  * Hollenbach et al. (1994) Hollenbach, D., Johnstone, D., Lizano, S., & Shu, F. 1994, ApJ, 428, 654, doi: 10.1086/174276
  * Hunter (2007) Hunter, J. D. 2007, Computing In Science & Engineering, 9, 90
  * Jayawardhana et al. (2003) Jayawardhana, R., Mohanty, S., & Basri, G. 2003, ApJ, 592, 282, doi: 10.1086/375573
  * Johnstone et al. (2021) Johnstone, C. P., Bartel, M., & Güdel, M. 2021, A&A, 649, A96, doi: 10.1051/0004-6361/202038407
  * Kastner et al. (2016) Kastner, J. H., Principe, D. A., Punzi, K., et al. 2016, AJ, 152, 3, doi: 10.3847/0004-6256/152/1/3
  * Kounkel & Covey (2019) Kounkel, M., & Covey, K. 2019, AJ, 158, 122, doi: 10.3847/1538-3881/ab339a
  * Kuchner et al. (2016) Kuchner, M. J., Silverberg, S. M., Bans, A. S., et al. 2016, ApJ, 830, 84, doi: 10.3847/0004-637X/830/2/84
  * Lee et al. (2020) Lee, J., Song, I., & Murphy, S. 2020, MNRAS, 494, 62, doi: 10.1093/mnras/staa689
  * Mann et al. (2015) Mann, R. K., Andrews, S. M., Eisner, J. A., et al. 2015, ApJ, 802, 77, doi: 10.1088/0004-637X/802/2/77
  * Mann et al. (2014) Mann, R. K., Di Francesco, J., Johnstone, D., et al. 2014, ApJ, 784, 82, doi: 10.1088/0004-637X/784/1/82
  * Murphy et al. (2018) Murphy, S. J., Mamajek, E. E., & Bell, C. P. M. 2018, MNRAS, 476, 3290, doi: 10.1093/mnras/sty471
  * Nakatani et al. (2018) Nakatani, R., Hosokawa, T., Yoshida, N., Nomura, H., & Kuiper, R. 2018, ApJ, 865, 75, doi: 10.3847/1538-4357/aad9fd
  * Nicholson et al. (2019) Nicholson, R. B., Parker, R. J., Church, R. P., et al. 2019, MNRAS, 485, 4893, doi: 10.1093/mnras/stz606
  * Owen et al. (2012) Owen, J. E., Clarke, C. J., & Ercolano, B. 2012, MNRAS, 422, 1880, doi: 10.1111/j.1365-2966.2011.20337.x
  * Owen et al. (2010) Owen, J. E., Ercolano, B., Clarke, C. J., & Alexander, R. D. 2010, MNRAS, 401, 1415, doi: 10.1111/j.1365-2966.2009.15771.x
  * Parker et al. (2021) Parker, R. J., Alcock, H. L., Nicholson, R. B., Panić, O., & Goodwin, S. P. 2021, ApJ, 913, 95, doi: 10.3847/1538-4357/abf4cc
  * Pascucci et al. (2020) Pascucci, I., Banzatti, A., Gorti, U., et al. 2020, ApJ, 903, 78, doi: 10.3847/1538-4357/abba3c
  * Pecaut & Mamajek (2013) Pecaut, M. J., & Mamajek, E. E. 2013, ApJS, 208, 9, doi: 10.1088/0067-0049/208/1/9
  * Pecaut & Mamajek (2016) —. 2016, MNRAS, 461, 794, doi: 10.1093/mnras/stw1300
  * Picogna et al. (2021) Picogna, G., Ercolano, B., & Espaillat, C. C. 2021, MNRAS, 508, 3611, doi: 10.1093/mnras/stab2883
  * Picogna et al. (2019) Picogna, G., Ercolano, B., Owen, J. E., & Weber, M. L. 2019, MNRAS, 487, 691, doi: 10.1093/mnras/stz1166
  * Preibisch et al. (2005) Preibisch, T., Kim, Y.-C., Favata, F., et al. 2005, ApJS, 160, 401, doi: 10.1086/432891
  * Rieke et al. (2005) Rieke, G. H., Su, K. Y. L., Stansberry, J. A., et al. 2005, ApJ, 620, 1010, doi: 10.1086/426937
  * Scally & Clarke (2001) Scally, A., & Clarke, C. 2001, MNRAS, 325, 449, doi: 10.1046/j.1365-8711.2001.04274.x
  * Schneider et al. (2019) Schneider, A. C., Shkolnik, E. L., Allers, K. N., et al. 2019, AJ, 157, 234, doi: 10.3847/1538-3881/ab1a26
  * Silverberg et al. (2021) Silverberg, S. M., Günther, H. M., Kim, J. S., Principe, D. A., & Wolk, S. J. 2021, AJ, 162, 279, doi: 10.3847/1538-3881/ac2cc0
  * Silverberg et al. (2016) Silverberg, S. M., Kuchner, M. J., Wisniewski, J. P., et al. 2016, ApJ, 830, L28, doi: 10.3847/2041-8205/830/2/L28
  * Silverberg et al. (2018) —. 2018, ApJ, 868, 43, doi: 10.3847/1538-4357/aae3e3
  * Silverberg et al. (2020) Silverberg, S. M., Wisniewski, J. P., Kuchner, M. J., et al. 2020, ApJ, 890, 106, doi: 10.3847/1538-4357/ab68e6
  * Skinner & Güdel (2013) Skinner, S. L., & Güdel, M. 2013, ApJ, 765, 3, doi: 10.1088/0004-637X/765/1/3
  * Störzer & Hollenbach (1999) Störzer, H., & Hollenbach, D. 1999, ApJ, 515, 669, doi: 10.1086/307055
  * Tanaka et al. (2013) Tanaka, K. E. I., Nakamoto, T., & Omukai, K. 2013, ApJ, 773, 155, doi: 10.1088/0004-637X/773/2/155
  * Telleschi et al. (2007) Telleschi, A., Güdel, M., Briggs, K. R., Audard, M., & Palla, F. 2007, A&A, 468, 425, doi: 10.1051/0004-6361:20066565
  * Theissen & West (2014) Theissen, C. A., & West, A. A. 2014, ApJ, 794, 146, doi: 10.1088/0004-637X/794/2/146
  * Tody (1993) Tody, D. 1993, in Astronomical Society of the Pacific Conference Series, Vol. 52, Astronomical Data Analysis Software and Systems II, ed. R. J. Hanisch, R. J. V. Brissenden, & J. Barnes, 173
  * van Terwisga et al. (2019) van Terwisga, S. E., Hacar, A., & van Dishoeck, E. F. 2019, A&A, 628, A85, doi: 10.1051/0004-6361/201935378
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261, doi: https://doi.org/10.1038/s41592-019-0686-2
  * Wahhaj et al. (2010) Wahhaj, Z., Cieza, L., Koerner, D. W., et al. 2010, ApJ, 724, 835, doi: 10.1088/0004-637X/724/2/835
  * Wang & Goodman (2017) Wang, L., & Goodman, J. 2017, ApJ, 847, 11, doi: 10.3847/1538-4357/aa8726
  * Wilhelm & Portegies Zwart (2022) Wilhelm, M. J. C., & Portegies Zwart, S. 2022, MNRAS, 509, 44, doi: 10.1093/mnras/stab2523
  * Williams & Cieza (2011) Williams, J. P., & Cieza, L. A. 2011, ARA&A, 49, 67, doi: 10.1146/annurev-astro-081710-102548
  * Winter et al. (2019) Winter, A. J., Clarke, C. J., & Rosotti, G. P. 2019, MNRAS, 485, 1489, doi: 10.1093/mnras/stz473
  * Wölfer et al. (2019) Wölfer, L., Picogna, G., Ercolano, B., & van Dishoeck, E. F. 2019, MNRAS, 490, 5596, doi: 10.1093/mnras/stz2939
  * Wyatt (2008) Wyatt, M. C. 2008, ARA&A, 46, 339, doi: 10.1146/annurev.astro.45.051806.110525
  * Yang et al. (2012) Yang, H., Herczeg, G. J., Linsky, J. L., et al. 2012, ApJ, 744, 121, doi: 10.1088/0004-637X/744/2/121

##  Appendix A Best-Fit Spectral Analysis Models

We describe our spectral fitting procedure in Section 2. The resulting best
fits from this analysis are shown in Figures 5-8, respectively. We discuss the
derived X-ray properties of our Peter Pan sources in Section 3.

##  Appendix B M star X-ray Control Sample

To gauge the extent to which the Peter Pan systems are potentially X-ray
under-luminous, we construct a control sample of objects to compare against.
For an appropriate comparison, we only consider objects with suspected
spectral types in a similar range as that identified for the Peter Pan sources
from Silverberg et al. (2020) (i.e. M4–M7). We also require high
($\gtrsim$95%) membership probability in the suspected parent moving groups of
the identified Peter Pan sources (i.e. Carina, Columba, and Tuc-Hor) as
determined by the BAYNAN algorithm (Gagné et al., 2018). These populations are
expected to be co-eval with derived ages of $\sim$45 Myr (Bell et al., 2015).

For this initial sample of 227 objects (J. Gagné, private communication), we
simultaneously query the XMM-Newton Serendipitous Source Catalog (4XMM-DR11)
as well as the Chandra Source Catalog (v.2.0) for XMM/Chandra pipeline-derived
X-ray fluxes reported by the HEASARC interface. We consider a coordinate
search radius of 5″ for each object, returning 13 matches. We carefully
inspect the X-ray light curves of these matches by eye when available,
excluding one source due to an obvious flare feature. To compute X-ray
luminosities, we use respective Gaia DR3 distances along with the pipeline
derived fluxes listed (EP 8 band for XMM-Newton matches, CSC broadband for
Chandra matches). We note these energy ranges (0.2-12 keV and 0.5-7.0 keV,
respectively) do not exactly line up with our reported 0.3-8 keV luminosities.
From our best fit spectral models, we compute our Peter Pan X-ray luminosities
in this larger energy band (0.2-12 keV) but find the overall difference is
negligible, given the X-ray flux strongly peaks at 1 keV. For this sub-sample
of objects, we also compute bolometric luminosities as indicated in Section 2.
We note we do not consider X-ray fluxes derived by ROSAT observations, given
its lower overall sensitivity biases detections towards the brightest, X-ray
luminous objects. We report the overall X-ray properties of this control
sample in Table 5.

Table 5: Analog M Star X-ray Sample Source | Moving GroupaaThe Tuc-Hor and Carina moving groups are estimated to be 45 Myr old (Bell et al., 2015). | Spectral Type | Lx (erg/s) | log Lx/Lbol
---|---|---|---|---
2MASS J01275875-6032243 | THA | M4 | 2.3e+28 | -3.33
2MASS J02025502-7259155 | THA | M5.5 | 2.87e+28 | -3.2
2MASS J02153328-5627175 | THA | M4.5 | 9.14e+27 | -3.28
2MASS J03152363-5342539 | THA | M5.2 | 1.9e+28 | -3.27
WOH G 312 | CAR | M4.6 | 1.96e+28 | -3.42
2MASS J08420090-7113216 | CAR | M4.7 | 1.44e+29 | -2.81
2MASS J09421385-5601361 | CAR | M4 | 2.34e+29 | -2.63
UCAC3 101-434911 | THA | M4 | 4.06e+28 | -3.26
2MASS J23143092-5405313 | THA | M5 | 1.12e+28 | -3.3
2MASS J23143092-5405313 | THA | M5 | 4.71e+27 | -3.68
2MASS J05360322-6555191 | THA | M6.5 | 3.27e+27 | -3.4
Figure 5: Chandra/ACIS X-ray spectra of J0949A (left) and J0949B (right)
observed on 2019-12-13. Blue points represent our binned X-ray count data. The
orange line represents our best-fit spectral model (see Section 2). The data
is fit without binning using the Cstat statistic, but is shown here binned to
5 counts/bin for display purposes (error bars are $\sqrt{N}$) to guide the
eye. Figure 6: Similar to Figure 5, Chandra/ACIS X-ray spectra of J0501
observed on 2019-12-22. Figure 7: Similar to Figure 5, Chandra/ACIS X-ray
spectra of J0446A (top row) and J0446B (bottom row) observed on 2020-05-20
(left), 2020-05-21 (middle) and 2020-06-02 (right). Figure 8: Similar to
Figure 5, Chandra/ACIS X-ray spectra of J0808 observed on 2021-07-06 (left)
and 2021-07-07 (right).


# Sharding-Based Proof-of-Stake Blockchain Protocols: Security Analysis

Abdelatif Hafid, Abdelhakim Senhaji Hafid, Adil Senhaji, A. Hafid and A. S.
Hafid are with Montreal Blockchain Laboratory (mbl), Department of Computer
Science and Operational Research, University of Montreal, Montreal, QC H3T
1J4, Canada.
E-mail<EMAIL_ADDRESS>A. Senhaji is a Managing Director of
technology at Mizuho Securities, USA.

###### Abstract

Blockchain technology has been gaining great interest from a variety of
sectors, including healthcare, supply chain and cryptocurrencies. However,
Blockchain suffers from its limited ability to scale (i.e. low throughput and
high latency). Several solutions have been appeared to tackle this issue. In
particular, sharding proved that it is one of the most promising solutions to
Blockchain scalability. Sharding can be divided into two major categories: (1)
Sharding-based Proof-of-Work (PoW) Blockchain protocols, and (2) Sharding-
based Proof-of-Stake (PoS) Blockchain protocols. The two categories achieve a
good performances (i.e. good throughput with a reasonable latency), but raise
security issues. This article attends that analyze the security of the second
category. More specifically, we compute the probability of committing a faulty
block and measure the security by computing the number of years to fail.
Finally, to show the effectiveness of the proposed model, we conduct a
numerical analysis and evaluate the results obtained.

###### Index Terms:

blockchain scalability, sharding, security analysis, Proof-of-Stake, practical
Byzantine fault tolerance

##  1 Introduction

With the raise of Bitcoin [1], Blockchain has attracted significant attention
and extensive research. More specifically, it has been adopted by numerous and
several applications, such as healthcare [2], supply chain [3], and Internet-
of-Things [4]. However, Blockchain’s capacity to scale is very limited [5].
For example, in the case of cryptocurrencies, Bitcoin [1] handles between 3-7
transactions per second (tps), which is very limited compared to traditional
payment systems (e.g. PayPal [6]). To deal with this issues, many solutions
appeared [5]. In particular, sharding is emerged as a promising solution [5].
Sharding consists of partitioning the network into shanks, called shards; all
shards work in parallel to enhance the performance of the network and then
mitigate the scalability issues. More specifically, each shard handles a sub-
set of transactions instead of the entire network handles all the
transactions. However, the security of sharding-based blockchain protocols is
emerging as a challenging issue. More specifically, in sharding-based
blockchain environment, it is easy for a malicious user to conquer and attack
a single shard compared to the whole network. This attack is well-known as a
shard takeovers attack [7].

Malicious nodes (e.g. Sybil nodes) are increasingly appear with the growing
and the spread of Blockchain technology. In contract, several consensus
mechanisms appeared to deal with these Byzantine nodes, including, PoW, PoS, ,
and practical Byzantine Fault Tolerance (pBFT).

pBFT is an algorithm that tolerates Byzantine fault [8]; this algorithm
belongs to the BFT class. Byzantine Fault Tolerance (BFT) is the feature of a
distributed network to reach consensus (agreement on the same value) even when
some of nodes in the network fail to respond or respond with incorrect
information. Leslie Lamport [9] proved that if we have $3m+1$ correctly
working processors, a consensus (agreement on same value) can be reached if at
most $m$ processors are faulty, which means strictly more than two-thirds of
the total number of processors should be honest.

PoS is an alternative consensus mechanism of PoW. It consists of selecting
validators in proportion to their number of coins. Validators are responsible
of adding new blocks to Blockchain.

By going through the literature [10, 5], we can classify sharding-based
Blockchain protocols into two major categories: sharding-based PoW and
sharding-based PoS Blockchain protocols.

Recently, Hafid et al. [11, 7, 12] proposed many approaches and models to
analyze the security of the first category. However, there is a lack of
methods in the literature to analyze the second category. In this paper, we
focus on the second category. Precisely, we provide a probabilistic model to
analyze the security of this kind of sharding-based Blockchain protocols by
computing the probability of committing a faulty block. Furthermore, based on
these probabilities, we calculate the number of years to fail for the purpose
of quantifying and measuring the security of the network.

The rest of the paper is structured as follows. Section 2 presents an overview
of sharding-based PoS Blockchain protocols and presents the proposed
probabilistic model. Section 3 presents numerical results and evaluates the
proposed model. In the end, we conclude this paper in Section 4.

##  2 Mathematical Model

The objective of this section is to propose a probabilistic model to analyze
the security of sharding-based PoS Blockchain protocols.

###  2.1 Abbreviations and Definitions

In what follows, we provide abbreviations and definitions required for the
rest of the paper. Table I shows the list of symbols and variables that are
used to describe the proposed model.

TABLE I: Notations & Symbols Notation |  Description
---|---
$N$ |  Number of users
$n$ |  Size of a shard’s committee
$n^{\prime}$ |  Size of the beacon’s committee
$H$ |  Number of honest validators in a shard
$M$ |  Number of malicious validators in a shard
$V$ |  Number of validators in a shard ($V=H+M$)
$\zeta$ |  Number of committees
$X$ |  Random variable that computes the number of malicious nodes in the committee of a shard
$H^{\prime}$ |  Number of honest validators in the beacon chain
$M^{\prime}$ |  Number of malicious validators in the beacon chain
$V^{\prime}$ |  Number of validators in the beacon chain ($V^{\prime}=H^{\prime}+M^{\prime}$)
$X^{\prime}$ |  Random variable that computes the number of malicious nodes in the committee of the beacon chain
$r$ |  Committee/beacon resiliency
$R$ |  Shard resiliency
$R^{\prime}$ |  Beacon resiliency
$\mathcal{P}$ |  Probability of conquering the protocol
$P$ |  Probability of a shard to commit a faulty block
$P^{{}^{\prime}}$ |  Probability of the beacon chain to commit a faulty block
$P^{{}^{\prime\prime}}$ |  Probability of all shards committing a faulty block
$p_{m}$ |  Percentage of malicious validators in a shard chain as well as in the beacon chain
$Y_{f}$ |  Years to fail

######  Definition 2.1 (Shard’s committee).

Shard’s committee is a subset of validators selected randomly from the set of
validators that are decided to stake for the shard chains.

######  Definition 2.2 (Beacon’s committee).

Beacon’s committee is a subset of validators selected randomly from the set of
validators that are decided to stake for the beacon chain.

###  2.2 Probabilistic Model

In this section, we present the proposed probabilistic model. Figure 1 shows a
network of sharding-based PoS Blockchain protocol. The network contains a
single beacon chain and $\zeta$ shard chains. Shard chains forge/produce
blocks in parallel. All shard chains are synchronized by beacon chain. More
specifically, each shard has its own committee, which is randomly assigned by
the beacon chain. And each shard committee processes transactions belonging to
it. When a shard block is created, the beacon committee verify the block, if
valid, it adds the block header to the beacon chain. Otherwise, it drops it
and sends the proof to other shards to vote for the purpose of slashing the
misbehaving shard committee. Furthermore, in each epoch (a period of time),
beacon chain shuffles committees to ensure the security. For Incognito [13],
when a new random number is generated, the beacon chain shuffles committees;
this means that, one epoch for Incognito is corresponding to generating a new
random number. This number is generated periodically in a round-robin fashion
[13], [14].

Shard 1Shard nBeacon Chain. . . .Honest NodeByzantine NodeUsers (nodes)
staketo be validators to the Beacon chainUsers (nodes) staketo be validators
to the shards chains Figure 1: A scheme of a sharding-based PoS and pBFT
Blockchain protocol

######  Lemma 2.1.

The probability of a shard’s committee members committing a faulty block ($P$)
can be described as follows:

$P(X\geq\frac{2n}{3})=\sum_{k=\frac{2n}{3}}^{n}\frac{\binom{M}{k}\binom{H}{n-k}}{\binom{V}{n}}$
(1)

######  Lemma 2.2.

The probability of at least $\frac{2}{3}$ of all shards committees committing
a faulty block ($P^{{}^{\prime}}$) can be computed as follows:

$\sum_{i=\frac{2\zeta}{3}}^{\zeta}\bigg{(}P(X\geq\frac{2n}{3})\bigg{)}^{i}=\sum_{i=\frac{2\zeta}{3}}^{\zeta}\sum_{\alpha=\frac{2n}{3}}^{n}\Bigg{(}\frac{\binom{M}{\alpha}\binom{H}{n-\alpha}}{\binom{V}{n}}\Bigg{)}^{i}$
(2)

###### Proof.

The minimum of all shards’ committees to commit a faulty block is
$\frac{2\zeta}{3}$, where $\zeta$ is the number of shards. The probability of
exactly $\frac{2\zeta}{3}$ all shards’ committees confirm/agree to add a
faulty block can be summarized and expressed as follows:

$P_{\frac{2\zeta}{3}}=\bigg{(}P(X\geq\frac{2n}{3})\bigg{)}^{\frac{2\zeta}{3}}$
(3)

And the probability to commit a faulty block by exactly $\frac{2\zeta}{3}+1$
of all shards’ committees can be described as follows.

$P_{\frac{2\zeta}{3}+1}=\bigg{(}P(X\geq\frac{2n}{3})\bigg{)}^{\frac{2\zeta}{3}+1}$
(4)

Similarly, the probability of exactly $\zeta$ of all shards’ committees (the
entire number of shards in this case) agree to add a faulty block can be
expressed as follows:

$P_{\zeta}=\bigg{(}P(X\geq\frac{2n}{3})\bigg{)}^{\zeta}$ (5)

A faulty block can be committed if $\frac{2\zeta}{3}$ or $\frac{2\zeta}{3}+1$
or $\frac{2\zeta}{3}+2$ or, $\cdots$, or $\zeta$ of all shards’ committees
agree to add this block. This can be mathematically computed by the sum over
all these probabilities and can be expressed and summarized as follows:

$P^{{}^{\prime}}=P_{\frac{2\zeta}{3}}+P_{\frac{2\zeta}{3}+1}+\cdots+P_{\zeta}$
(6)

∎

######  Lemma 2.3.

The probability of the beacon’s committee members committing a faulty block
($P^{{}^{\prime\prime}}$) can be expressed as follows:

$P(X^{\prime}\geq\frac{2n^{\prime}}{3})=\sum_{j=\frac{2n^{\prime}}{3}}^{n^{\prime}}\frac{\binom{M^{\prime}}{j}\binom{H^{\prime}}{n^{\prime}-j}}{\binom{V^{\prime}}{n^{\prime}}}$
(7)

Proofs of Lemma 2.1 and 2.2 are a direct results from the cumulative
hypergeometric distribution [12, 7].

######  Theorem 2.1 (Committing a Faulty Block).

The probability of committing a faulty block for a given shard can be
expressed as follows:

$\mathcal{P}=\sum_{k=\frac{2n}{3}}^{n}\sum_{i=\frac{2\zeta}{3}}^{\zeta}\sum_{\alpha=\frac{2n}{3}}^{n}\sum_{j=\frac{2n^{\prime}}{3}}^{n^{\prime}}\frac{\binom{M}{k}\binom{H}{n-k}\binom{M}{\alpha}^{i}\binom{H}{n-\alpha}^{i}\binom{M^{\prime}}{j}\binom{H^{\prime}}{n^{\prime}-j}}{\binom{V}{n}\binom{H}{n-\alpha}^{i}\binom{V^{\prime}}{n^{\prime}}}$
(8)

###### Proof.

To commit a faulty block, it must be confirmed/verified by at least
$\frac{2}{3}$ of a shard’s committee members and by at least $\frac{2}{3}$ of
beacon’s committee members, and finally by at least $\frac{2}{3}$ of all
shards’ committees. This can be identified mathematically by the product over
the three probabilities (the calculated probabilities in Lemmas 2.1, 2.2, and
2.3). ∎

##  3 Results & Evaluation

In Figures 2, 4, and 3, we assume a network with N = 2000 users/nodes, V =
200, V’= 400, r = 0.333, R = 0.1, R = 0.2, and R = 0.3. The remaining
parameters will be shown in Figures.

Figure 2: Log-scale plot of the probability of a shard to commit a faulty
block ($P$) versus the size of the committee (n)

Figure 2 shows the probability of a shard to commit a faulty block versus the
size of the committee in a network of 2000 nodes. We observe that the
probability $P$ decreases when the size of the committee increases. More
specifically, we observe that the probability corresponding to $R=0.1$ (i.e.
10% of malicious nodes in each shards) decreases rapidly compared to those of
$R=0.15$ and $R=0.2$; this can be explained by the small percentage of
malicious nodes. In other words, as the percentage of malicious nodes gets
smaller the probability decreases and vice versa.

Figure 3: Log-scale plot of the probability of all shards committing a faulty
block ($P^{{}^{\prime}}$) versus the size of the committee (n)

Figure 4 shows the probability of all shards committing a faulty block versus
the size of the committee. We observe that the probability $P^{{}^{\prime}}$
decreases when the size of the committee increases.

Figure 4: Log-scale plot of the probability of the beacon chain to commit a
faulty block ($P^{{}^{\prime\prime}}$) versus the size of the committee (n)

Figure 3 shows the probability of the beacon chain to commit a faulty block
versus the size of the committee. We observe that the probability
$P^{{}^{\prime\prime}}$ decreases when the size of the committee increases.
More specifically, we observe that the probability corresponding to $R=0.1$
(i.e. 10% of malicious nodes in each shards) decreases sharply compared to
those of $R=0.15$ and $R=0.2$.

TABLE II: Probability of conquering the protocol $p_{m}$ | 10 % | 15 % | 20 % | 30 %
---|---|---|---|---
$\mathcal{P}^{\mathrm{a}}$ | 3.63E-66 | 2.10E-34 | 1.58E-18 | 1.70E-04
$Y_{f}^{\mathrm{a}}$ | 7.56E+62 | 1.30E+31 | 1.74E+17 | 16.12
$\mathcal{P}^{\mathrm{b}}$ | 0.0 | 5.14E-80 | 2.01E-41 | 5.30E-07
$Y_{f}^{\mathrm{b}}$ | inf | 5.33E+76 | 1.36E+38 | 5171.32
aScenario 1; bScenario 2. |

In Figure II, we assume two scenarios to show the effectiveness and the
feasibility of the proposed model: Scenario 1 proposes a network with
$N=2000$, $\zeta=8$, $V=200$, $V^{\prime}=400$, and $r=r^{\prime}=0.333$
whereas Scenario 2 proposes a network with $N=4000$, $\zeta=8$, $V=400$,
$V^{\prime}=800$, and $r=r^{\prime}=0.333$. It is noteworthy that the proposed
model can be adopted to any scenario.

Table II shows the probability of conquering the chain (i.e. the probability
of committing a faulty block; it is calculated based on Theorem 2.1) for
different percentage of malicious nodes in the shards as well as in the beacon
chain. Moreover, Table II shows the number of years to fail corresponding to
these probabilities. More specifically, we observe that as the percentage of
malicious nodes increases the number of years to fail decreases.

We observe that the probability of conquering the chain is extremely low even
with 20% of malicious nodes in each shard as well as in the beacon chain. This
achieves a good security, which is about $1.74E+17$ years to fail.

To sum up, we conclude that by adjusting the size of the shard’s committee as
well as the size of the beacon’s committee, we could support/protect the
sharded Blockchain systems against malicious nodes (e.g. Sybil nodes).

##  4 Conclusion

This article addresses the security of sharding-based blockchain protocols
that are based on PoS and pBFT consensuses. In particular, it provides a
probabilistic model to compute the probability of committing a faulty block.
Based on this probability, we compute the number of years to fail.
Furthermore, this article depicts that we can control the number of years to
fail by adjusting the size of the shard’s committee as well as the size of the
beacon’s committee. Future works focus on computing the failure probability
across-shard transaction.

## References

  * [1] S. Nakamoto, “Bitcoin whitepaper,” _URL: https://bitcoin.org/bitcoin.pdf_ , 2008.
  * [2] M. H. Kassab, J. DeFranco, T. Malas, P. Laplante, g. destefanis, and V. V. Graciano Neto, “Exploring research in blockchain for healthcare and a roadmap for the future,” _IEEE Transactions on Emerging Topics in Computing_ , pp. 1–1, 2019.
  * [3] M. Du, Q. Chen, J. Xiao, H. Yang, and X. Ma, “Supply chain finance innovation using blockchain,” _IEEE Transactions on Engineering Management_ , vol. 67, no. 4, pp. 1045–1058, 2020.
  * [4] M. S. Ali, M. Vecchio, M. Pincheira, K. Dolui, F. Antonelli, and M. H. Rehmani, “Applications of blockchains in the internet of things: A comprehensive survey,” _IEEE Communications Surveys Tutorials_ , vol. 21, no. 2, pp. 1676–1717, 2019.
  * [5] A. Hafid, A. S. Hafid, and M. Samih, “Scaling blockchains: A comprehensive survey,” _IEEE Access_ , vol. 8, no. 19800223, pp. 125 244–125 262, 2020\.
  * [6] PayPal, 2021. [Online]. Available: https://www.paypal.com/fr/webapps/mpp/home
  * [7] A. Hafid, A. S. Hafid, and M. Samih, “A tractable probabilistic approach to analyze sybil attacks in sharding-based blockchain protocols,” _arXiv_ , 2020\. [Online]. Available: https://arxiv.org/abs/2104.07215
  * [8] M. Castro, B. Liskov _et al._ , “Practical byzantine fault tolerance,” in _OSDI_ , vol. 99, no. 1999, 1999, pp. 173–186.
  * [9] L. Lamport, R. Shostak, and M. Pease, “The byzantine generals problem,” in _Concurrency: the Works of Leslie Lamport_ , 2019, pp. 203–226.
  * [10] G. Wang, Z. J. Shi, M. Nixon, and S. Han, “Sok: Sharding on blockchain,” in _Proceedings of the 1st ACM Conference on Advances in Financial Technologies_ , 2019, pp. 41–61.
  * [11] A. Hafid, A. S. Hafid, and M. Samih, “New mathematical model to analyze security of sharding-based blockchain protocols,” _IEEE Access_ , vol. 7, no. 19232899, pp. 185 447–185 457, 2019.
  * [12] ——, “A novel methodology-based joint hypergeometric distribution to analyze the security of sharded blockchains,” _IEEE Access_ , vol. 8, no. 20000968, pp. 179 389–179 399, 2020.
  * [13] Incognito, 2021. [Online]. Available: https://we.incognito.org/t/incognito-whitepaper-incognito-mode-for-cryptonetworks/168
  * [14] R. V. Rasmussen and M. A. Trick, “Round robin scheduling–a survey,” _European Journal of Operational Research_ , vol. 188, no. 3, pp. 617–636, 2008.


# On the stability of equilibria of the physiologically-informed dynamic
causal model

Sayan Nag
<EMAIL_ADDRESS>
###### Abstract

Experimental manipulations perturb the neuronal activity. This phenomenon is
manifested in the fMRI response. Dynamic causal model and its variants can
model these neuronal responses along with the BOLD responses [1, 2, 3, 4, 5] .
Physiologically-informed DCM (P-DCM) [5] gives state-of-the-art results in
this aspect. But, P-DCM has more parameters compared to the standard DCM model
and the stability of this particular model is still unexplored. In this work,
we will try to explore the stability of the P-DCM model and find the ranges of
the model parameters which make it stable.

##  1 Introduction

Experimental perturbations lead to the fluctuations in the neuronal activity
which are known as neuronal responses. These neuronal responses drive the
vasoactive signals which in turn drive the blood flow responses. The changes
in the blood flow responses lead to changes in the blood volume and
deoxyhemoglobin content which together give rise to BOLD signals.

To model this entire flow which occurred in the brain due to some activity and
find the causal interactions between the different inter-connected brain
regions biophysical models like Dynamic Causal Models (DCMs) have been
generated. DCM explains the causality between different regions. The
connection between different regions are given in the form of a probabilistic
graphical model where the nodes represent the brain regions. The effective
connectivity between these regions are parameterized in terms of the coupling
parameters. In other words, the effective connectivity between two regions or
nodes is expressed in terms of the conditional dependency between those two
regions.

The physiologically-informed DCM (P-DCM) is a physiologically realistic state-
of-the-art DCM variant which gives better model evidence when compared to
other DCM approaches. It introduces a factor in the connection between the
excitatory and inhibitory neuronal population. Along with this, there are some
modifications in the neurovascular coupling where a feed-forward model has
been introduced instead of feedback-based model. It also introduces a dynamic
viscoelastic effect in the balloon model. Altogether, it introduces more
parameters compared to the vanilla DCM, renders dynamic stability to the
system and makes the model more realistic from physiological point of view.

Equilibrium points are not always stable. It is important to analyze the
stability of the equilibrium points they play different roles in determining
the dynamics of a system. In this work, considering each of the model
equations we will look into the parameters and their corresponding ranges
which make the pDCM model dynamically stable near respective equilibrium
points.

##  2 P-DCM Model

The P-DCM model consists of 3 sub-models, namely, the neuronal model, the
neurovascular coupling, the hemodynamic and balloon models and the BOLD model.
The BOLD model consists a set of equations which is empirically determined. We
will focus on the neuronal model, the neurovascular coupling and the
hemodynamic and balloon models which were the main contributions of the P-DCM
paper [5].

###  2.1 Neuronal model

The neuronal model consists of the excitatory and inhibitory states. The
dynamics are dictated by the following state-space differential equations:

$\frac{dx_{E}(t)}{dt}=-\sigma x_{E}(t)-\mu x_{I}(t)+cu(t)$ (1)
$\frac{dx_{I}(t)}{dt}=\lambda x_{E}(t)-\lambda x_{I}(t)$ (2)

Here, $\sigma$ represents the self-intrinsic connectivity of the excitatory
neuron. $\lambda$ is the inhibitory gain factor. The excitatory-inhibitory
mutual connection is given by $\mu$.

###  2.2 Neurovascular coupling (NVC)

Neurovascular coupling refers to the changes in CBF and CBV due to change in
the neuronal activity. A feedforward based NVC model has been chosen in P-DCM
which is given by the following set of differential equations:

$\frac{da(t)}{dt}=-\rho a(t)+x_{E}(t)$ (3) $\frac{df(t)}{dt}=\phi
a(t)-\chi(f(t)-1)$ (4)

Here, a(t) is the vasoactive signal. Using the above equations, it transforms
the neuronal response $x_{E}(t)$ to the blood flow response which is given by
f(t), $\rho$, $\phi$, $\chi$ respectively represents the decay of vasoactive
signal, the gain of vasoactive signal and the decay of blood inflow signal.

###  2.3 Hemodynamic and Balloon models

The hemodynamic model dictates the dynamics of the hemodynamic variables
called blood volume v(t) and the deoxyhemoglobin content q(t) with respect to
the blood inflow f(t) and the blood outflow $f_{out}(v)$.

$\frac{dv(t)}{dt}=\frac{1}{t_{MTT}}[f(t)-f_{out}(v,t)]$ (5)
$\frac{df(t)}{dt}=\frac{1}{t_{MTT}}[f(t)\frac{E(f)}{E_{0}}-f_{out}(v,t)\frac{q(t)}{v(t)}]$
(6)

Here, $t_{MTT}$ is the mean transit time the blood takes to pass the veins.
E(f) is the oxygen extraction fraction when the blood inflow is f(t) and
$E_{0}$ is the net oxygen extraction at rest. E(f) is given as:

$E(f)=1-(1-E_{0})^{1/f}$ (7)

Considering the visco-elastic effect the balloon model is given as:

$f_{out}(v,t)=\frac{1}{\tau+t_{MTT}}(t_{MTT}v(t)^{1/\alpha}+\tau f(t))$ (8)

Here, $\tau$ is the visco-elastic time constant and $\alpha$ is the Grubb’s
exponent.

##  3 Analysis of P-DCM Model

Linear Stability Analysis refers to the stability analysis of a system after
the linearization of the differential equations at respective equilibrium
points are done. Equilibrium refers to the state of a system which does not
change. Equilibrium point(s) can be estimated by setting the derivative(s) of
the differential equation(s) (describing the system) to zero.

A Jacobian matrix is formed by the first order partial derivatives of a system
of differential equations. The eigen values of a Jacobian matrix can be
represented in the form of complex numbers. The complex part of the eigenvalue
contributes an oscillatory component, but, the real part determines the
stability of the equilibrium point.

An equilibrium point of the differential equation defining a system is
considered to be: (a) stable if all the eigenvalues of the Jacobian matrix
evaluated at that equilibrium point have negative real parts and (b) unstable
if at least one of the eigenvalues of the Jacobian matrix evaluated at that
equilibrium point has a positive real part.

In this section we will show the detailed analysis of P-DCM model, especially
we will try to investigate the values of the model parameters which make the
model dynamically stable near the respective equilibrium points.

###  3.1 Neuronal model

Equations (1) and (2) give the neuronal dynamics of the P-DCM Model. Now,
making both of them equal to 0, we will get the equilibrium values of the
system. The equilibrium values of $x_{E}$ and $x_{I}$ are the same and equals
to $cu/(\sigma+\mu)$. Now, we will try to find the nature of stability of the
equilibrium point by computing the Jacobian matrix. The Jacobian matrix is
represented as:

$J_{NeuronalModel}=\begin{bmatrix}-\sigma&-\mu\\\
\lambda&-\lambda\end{bmatrix}$ (9)

The characteristic equation in terms of eigen values (s) is given as:

$s^{2}+(\sigma+\lambda)s+(\mu\lambda+\sigma\lambda)=0$ (10)

$s_{1}$ and $s_{2}$ are the roots of the above equation given as:

$s_{1}=-(\sigma+\lambda)+\sqrt{(\sigma+\lambda)^{2}-4(\mu\lambda+\sigma\lambda)}$
(11)

$s_{2}=-(\sigma+\lambda)-\sqrt{(\sigma+\lambda)^{2}-4(\mu\lambda+\sigma\lambda)}$
(12)

So, for the equilibrium point to be a stable one (and hence the system to be
stable), 1\. $\sigma+\lambda$ has to be more than 0. 2\.
$\sqrt{(\sigma+\lambda)^{2}-4(\mu\lambda+\sigma\lambda)}$ can be $\leq 0$, 3\.
$\sigma+\lambda>\mod{\sqrt{(\sigma+\lambda)^{2}-4(\mu\lambda+\sigma\lambda)}}$,
if $\sqrt{(\sigma+\lambda)^{2}-4(\mu\lambda+\sigma\lambda)}>0$

By definition each of $\sigma$ and $\lambda$ is positive. So the first
condition is satisfied. Considering $\lambda$ > 0,for the third condition, we
see that $\mu+\sigma$ has to be greater than 0. By, definition $\mu$ is also
positive, so third condition holds true. For the second condition we get:

$(\sigma-\lambda)^{2}\leq 4\mu\lambda$
$\lambda-2\sqrt{\mu\lambda}\leq\sigma\leq\lambda+2\sqrt{\mu\lambda}$

The given ranges for $\mu$ and $\lambda$ in the P-DCM paper [5] are:

$0<\mu<1.5,$ $0<\lambda<0.3$

Now, using these ranges for $\mu$ and $\lambda$ in the above expression
showing the permissible range of $\sigma$ and also considering $\sigma$ to be
positive (by definition and also to ensure stability), we find:

$0<\sigma<1.642$

In the P-DCM paper, $\sigma$ was made to vary within 0 and 1.5. So as already
computed above, the mentioned range of sigma (though it can be extended to
1.64) makes the equilibrium a stable focus and thus makes the system
asymptotically stable.

###  3.2 Neurovascular Coupling (NVC)

Equations (3) and (4) give the NVC dynamics of the P-DCM Model. Now, making
both of them equal to 0, we will get the equilibrium values of the system. The
equilibrium values of $a$ and $f$ are respectively $\frac{x_{E}^{*}}{\rho}$
and $\frac{\phi x_{E}^{*}}{\rho\chi}+1$. Here, $x_{E}^{*}$ is the steady-state
value of $x_{E}$. Now, we will try to find the nature of stability of the
equilibrium point by computing the Jacobian matrix. The Jacobian matrix is
represented as:

$J_{NVC}=\begin{bmatrix}-\rho&0\\\ \phi&-\chi\end{bmatrix}$ (13)

The characteristic equation in terms of eigen values (s) is given as:

$(s+\chi)(s+\rho)=0$ (14)

$s_{1}$ and $s_{2}$ are the roots of the above equation given as:

$s_{1}=-\chi$ (15)

$s_{2}=-\rho$ (16)

So, for the equilibrium point to be a stable one (and hence the system to be
stable), both $\chi$ and $\rho$ have to be positive. Based on the values
mentioned in the P-DCM paper ($\chi=0.6$ and $\rho=0.6$), the equilibrium is a
stable focus.

###  3.3 Hemodynamic and Balloon models

Putting the values of $E(f)$ and $f_{out}(v,t)$ in equations (5) and (6) and
making both of them equal to 0, we will get the equilibrium values of the
variables of the system. The equilibrium values of $v$ and $q$ are
respectively $v_{eq}$ and $q_{eq}$.

$v_{eq}=(f^{*})^{\alpha}$ (17)

$q_{eq}=\frac{v_{eq}f^{*}E^{*}}{f_{out}^{*}E_{0}}$ (18)

Here, $f^{*}$, $f_{out}^{*}$ and $E^{*}$ are respectively the steady-state
values of $f$, $f_{out}$ and $E$.

Now, we will try to find the nature of stability of the equilibrium point by
computing the Jacobian matrix. The Jacobian matrix is represented as:

$J_{H}=\begin{bmatrix}-\frac{v_{eq}^{1/\alpha-1}}{\alpha(t_{MTT}+\tau)}&0\\\
(\frac{1}{\alpha}-1)\frac{q_{eq}v_{eq}^{(1/\alpha-2)}}{t_{MTT}+\tau}+\frac{\tau
q_{eq}v_{eq}^{-2}}{t_{MTT}}&-(\frac{v_{eq}^{(1/\alpha-1)}}{t_{MTT}+\tau}+\frac{\tau
f^{*}}{t_{MTT}v_{eq}})\end{bmatrix}$ (19)

So, the eigen values $s_{1}$ and $s_{2}$ are given as:

$s_{1}=-\frac{v_{eq}^{1/\alpha-1}}{\alpha(t_{MTT}+\tau)}$ (20)

$s_{2}=-(\frac{v_{eq}^{(1/\alpha-1)}}{t_{MTT}+\tau}+\frac{\tau
f^{*}}{t_{MTT}v_{eq}})$ (21)

So, for the equilibrium point to be a stable one (and hence the system to be
stable), both $s_{1}$ and $s_{2}$ have to be negative. Based on the values
mentioned in the P-DCM paper ($1<t_{MTT}<5$ and $0<\tau<30$), the equilibrium
is a stable focus.

##  4 Conclusions

Linear stability analysis near the equilibrium points for the system of
differential equations defining the P-DCM model gives an idea about the ranges
of the parameters which make the model dynamically stable. It has been found
that the ranges of the model parameters defined in the actual P-DCM paper [5]
confer dynamical stability on the equilibrium points of the model differential
equations.

## References

  * [1] Karl J Friston, Lee Harrison, and Will Penny.  Dynamic causal modelling.  Neuroimage, 19(4):1273–1302, 2003.
  * [2] Klaas Enno Stephan, Lars Kasper, Lee M Harrison, Jean Daunizeau, Hanneke EM den Ouden, Michael Breakspear, and Karl J Friston.  Nonlinear dynamic causal models for fmri.  Neuroimage, 42(2):649–662, 2008.
  * [3] Rosalyn J Moran, Klaas E Stephan, T Seidenbecher, H-C Pape, Raymond J Dolan, and Karl J Friston.  Dynamic causal models of steady-state responses.  Neuroimage, 44(3):796–811, 2009.
  * [4] Jean Daunizeau, Karl J Friston, and Stefan J Kiebel.  Variational bayesian identification and prediction of stochastic nonlinear dynamic causal models.  Physica D: nonlinear phenomena, 238(21):2089–2118, 2009.
  * [5] Martin Havlicek, Alard Roebroeck, Karl Friston, Anna Gardumi, Dimo Ivanov, and Kamil Uludag.  Physiologically informed dynamic causal modeling of fmri data.  Neuroimage, 122:355–372, 2015.