11institutetext: Astrophysical Big Bang Laboratory, RIKEN, Saitama 351-0198,
Japan
email<EMAIL_ADDRESS>22institutetext: Interdisciplinary
Theoretical & Mathematical Science Program (iTHEMS), RIKEN, Saitama 351-0198,
Japan 33institutetext: Theoretical Astrophysics, IAAT, University of Tübingen,
72076 Tübingen, Germany 44institutetext: Department of Physics, Aristotle
University of Thessaloniki, Thessaloniki 54124, Greece

# Neutron star mass-radius constraints using the high-frequency QPOs of GRB
200415A

H. Sotani 1122 K. D. Kokkotas 3344 and N. Stergioulas 44

###### Abstract

Context. Quasi-periodic oscillations (QPOs) observed in a giant flare of a
strongly magnetized neutron star (magnetar), are carrying crucial information
for extracting the neutron star properties.

Aims. The aim of the study is to constrain the mass and radius of the neutron
star model for GRB 200415A, by identifying the observed QPOs with the crustal
torsional oscillations together with the experimental constraints on the
nuclear matter properties.

Methods. The frequencies of the crustal torsional oscillations are determined
by solving the eigenvalue problem with the Cowling approximation, assuming a
magnetic field of about $10^{15}$G.

Results. We find that the observed QPOs can be identified with several
overtones of crustal oscillations, for carefully selected combinations of the
nuclear saturation parameters. Thus, we can inversely constrain the neutron
star mass and radius for GRB 200415A by comparing them to the values of
nuclear saturation parameters obtained from terrestrial experiments.

Conclusions. We impose further constraints on the neutron star mass and radius
while the candidate neutron star models are consistent with the constraints
obtained from other available astronomical and experimental observations.

###### Key Words.:

asteroseismology – stars: neutron – stars: oscillations

##  1 Introduction

Neutron stars are remnants of core-collapse supernovae, in which extreme
states of matter are realized. The density inside the star easily exceeds the
standard nuclear density ($\rho_{0}$), while the gravitational and magnetic
fields inside/around the neutron star are among the strongest observed
anywhere in the Universe (Haensel et al. 2006). These extreme matter
conditions cannot be reproduced on Earth, and thus the observations of neutron
stars and the associated phenomena remain as the only way to extract
information about the state of matter at the most extreme densities. An
example of this “indirect” study of nuclear matter under extreme conditions is
the discovery of $2M_{\odot}$ neutron stars (Demorest et al. 2010; Antoniadis
et al. 2013; Cromartie et al. 2020). Based on these observations, the
possibility of so-called soft equations of state (EOS) for neutron star matter
has been practically abandoned (at least within the framework of Einstein’s
general relativity theory). In particular, many of the EOSs containing
hyperons, which generally have a soft core, have been excluded, leading to the
so-called hyperon puzzle (but see also e.g. Sun et al. (2023)). The detection
and analysis of the gravitational wave signal during the binary neutron star
merger GW170817 (Abbott et al. 2017, 2019), provided information on the tidal
deformability of the merging neutron stars, which can be translated to
constraints on the neutron star radius, see Chatziioannou (2020); Dietrich,
Hinderer & Samajdar (2021) for recent reviews. Furthermore, pulsar
observations provide useful constraints on the neutron star properties. For
example, the light bending due to the strong gravitational field induced by
the compact object, a purely relativistic effect, has been proven extremely
useful in constraining the neutron star properties. More specifically, since
the pulsar light curve depends on stellar compactness, one can constrain the
neutron star mass and radius by carefully analyzing the observed light curve
(e.g., (Pechenick et al. 1983; Leahy & Li 1995; Poutanen & Gierlinski 2003;
Psaltis & Özel 2014; Sotani & Miyamoto 2018)). Recently, the Neutron star
Interior Composition Explorer (NICER) on the International Space Station via
x-ray observations of neutron stars provided constraints on the parameters of
two neutron stars, PSR J0030+0451 (Riley et al. 2019; Miller et al. 2019) and
PSR J0740+6620 (Riley et al. 2021; Miller et al. 2021). Moreover, for the
first time, at the pre-merger stage preceding the GRB211211, quasi-periodic
oscillations (QPOs) at $22^{+3}_{-2}$ Hz and $51^{+2}_{-2}$Hz were detected
Xiao et al. (2022), suggesting that these features of the precursor signal may
have resulted from the resonant shattering (due to tidal interactions) of one
of the star’s crust prior to coalescence leading to the excitation of crustal
oscillations (Suvorov et al. 2022)111Note, however, that it is challenging to
reconcile the relatively short lifetime of a magnetar with a long inspiral
time of a binary neutron star system driven only by gravitational wave
radiation.

The modeling of the stellar oscillation spectrum provides a unique way for
studying the internal structure of stars. Especially for neutron stars, the
oscillations can be associated with the spectrum of the emitted gravitational
waves, and in this way, one can reveal details of their internal structure
(gravitational wave asteroseismology) (Andersson & Kokkotas 1996, 1998; Sotani
et al. 2004, 2011; Passamonti & Andersson 2012; Doneva et al. 2013; Krüger et
al. 2021). Moreover, gravitational wave asteroseismology has been used for
analyzing the gravitational signals in numerical simulations of the formation
of proto-neutron stars in core-collapse supernovae, e.g., Morozova et al.
(2018); Torres-Forné et al. (2019); Sotani et al. (2021).

In addition, the analysis of the afterglow spectra in magnetar giant flares
(SGR 1900+14 and SGR 1806-20) demonstrated a rich spectrum, the so-called
quasi-periodic oscillations (QPOs) (Strohmayer & Watts 2005; Israel et al.
2005; Strohmayer & Watts 2006). The analysis suggested that the QPOs were due
to magnetoelastic oscillation modes (Levin 2007; Sotani et al. 2007; Steiner &
Watts 2009; Gabler et al. 2011; Gearheart et al. 2011; Sotani et al. 2012;
Colaiuda & Kokkotas 2012; Gabler et al. 2012; Sotani et al. 2013; Gabler et
al. 2013, 2016, 2018) which provided another opportunity for constraining the
parameters of neutron stars.

Recently, another giant flare, GRB 200415A, was detected in the direction of
the NGC 253 galaxy by the Atmosphere-Space Interactions Monitor (ASIM) on the
International Space Station on the 15th of April 2020, where several high-
frequency QPOs, with varying significance, have been found at 836, 1444, 2132,
and 4250 Hz (Castro-Tirado et al. 2021). Unlike the QPOs observed in the
previous giant flares, the QPOs in GRB 200415A were confined only at high
frequencies, mainly due to the shortness of the observational interval.

Castro-Tirado et al. (2021) discuss two alternative explanations (that have
previously been proposed in the literature) for the occurrence of the high-
frequency QPOs in GRB 200415A. One possible explanation is that the high-
frequency QPOs are due to Alfvén waves traveling back and forth between the
footpoints of magnetic field lines, relatively close to the magnetar surface.
This may be the result of a reconnection event, following the development of
instability in the magnetosphere (Mahlmann et al. 2019). The second possible
explanation is the excitation of crustal oscillation, which is a viable
scenario for the QPOs reported for SGR 1900+14 and SGR 1806-20, as discussed
above. Due to the very short lifetime of the observed QPOs, Castro-Tirado et
al. (2021) slightly favors the first explanation, without excluding the
second.

Here, without attempting to favor one over the other of the two explanations,
we examine the consequences of the second one. That is, we try to answer the
question: If the high-frequency QPOs in GRB 200415A are due to crustal
oscillation, what constraints can be placed on the neutron star EOS?

Regarding the nature of the oscillations, one may attempt to associate some of
the observed high-frequency QPOs with polar-type oscillations of neutron
stars, such as the fundamental ($f$-) and pressure ($p$-) modes. Still, the
excitation of such global modes during giant flares is questionable, as they
require large density variations. It is thus more natural to associate these
oscillations with magnetoelastic or pure crustal torsional oscillations,
depending on the magnetic field strength.

The frequency and damping time of torsional oscillations in magnetars depends
on several factors, such as the equation of state, the properties of the
crust, the superfluidity of the core and the topology and strength of the
magnetic field, see Gabler et al. (2018) and references therein. Several
studies have shown that there are three distinct regions, in terms of the
effect of the magnetic field strength: a) At moderate magnetic field strengths
$\ll 10^{15}$G, the crustal oscillation frequencies are unaffected by the
magnetic field; b) For somewhat stronger fields of
$\,\,\raisebox{-3.44444pt}{$\stackrel{{\scriptstyle\textstyle<}}{{\sim}}$
}10^{15}$G, the frequencies are weakly affected, but the crustal oscillations
can be absorbed on a short timescale by a continuum of magneto-elastic
oscillations in the core222Notice that if the core is superfluid, then at
$\sim 10^{15}$G the Alfvén oscillations in the core can resonate with crustal
oscillations Gabler et al. (2013); Passamonti & Lander (2014); Gabler et al.
(2018). \- of particular relevance to the 3.5ms duration of the high-frequency
QPOs observed in GRB 200415A is the example of a timescale of $\sim 5$ms for
the damping of a crustal $n=1$ shear mode in a model with magnetic field
strength of $\sim 10^{15}$G in Gabler et al. (2018); c) At strengths
$\,\raisebox{-3.44444pt}{$\stackrel{{\scriptstyle\textstyle>}}{{\sim}}$
}10^{15}$G, the dominant character of the crustal oscillations becomes
magneto-elastic and gradually their frequency becomes strongly dependent on
the strength of the magnetic field.

Here, we will work under the assumption that the strength of the magnetic
field in GRB 200415A is in the second region, i.e. that it is
$\,\,\raisebox{-3.44444pt}{$\stackrel{{\scriptstyle\textstyle<}}{{\sim}}$
}10^{15}$G, so that it has a short damping time, but the frequencies are still
close to the values of pure crustal oscillations. This assumption is discussed
in more detail in Appendix A, where we show that for high-frequency torsional
oscillations, the shift of the crustal mode frequency caused by a magnetic
field with strength
$\,\,\raisebox{-3.44444pt}{$\stackrel{{\scriptstyle\textstyle<}}{{\sim}}$
}10^{15}$G is less than the observational uncertainty of the QPOs in GRB
200415A, as reported in Castro-Tirado et al. (2021) ($\sim 10\%$). In the
remainder of the paper, we will thus use frequencies of pure crustal torsional
oscillations (in the limit of no magnetic field).

In Sec. 2 we discuss the properties of the crust and their relation to
torsional oscillations. In Sec. 3 we try to identify the observed QPOs with
the higher overtones of torsional oscillations to arrive at EOS constraints.
The discussion of the results is resented in Sec. 4.

##  2 Crust equilibrium and torsional oscillations

Matter in the neutron star crust forms a Coulomb lattice, which behaves as a
solid (or liquid crystal). As a result, under perturbations the torsional
oscillations are favored. Since the matter in the core (inner or outer)
behaves as a nearly uniform density fluid the torsional oscillations are
confined only in the crust. Thus for a given EOS describing the core and the
crust and for a given central density, one can construct a neutron star model
by integrating the Tolman-Oppenheimer-Volkoff (TOV) equations starting with
appropriate initial conditions at the center. Here, we take a different
approach, specifying only the crust EOS and integrating the TOV equations
inwards (starting at the star’s surface), for various sets of mass $M$ and
radius $R$ (see Sotani et al. (2013) for details).

For the present study, we adopt the same phenomenological family of crust EOS
constructed by Oyamatsu & Iida (2003, 2007), the so-called OI-EOSs. For this
family of EOS the bulk energy per nucleon, for zero-temperature uniform
nuclear matter, is expressed in the vicinity of the saturation density for the
symmetric nuclear matter, $n_{0}$, as a function of the baryon number density,
$n_{\rm b}$, and the asymmetry parameter, $\alpha$, as

$w=w_{0}+\frac{K_{0}}{18n_{0}^{2}}(n_{\rm
b}-n_{0})^{2}+\left[S_{0}+\frac{L}{3n_{0}}(n_{\rm b}-n_{0})\right]\alpha^{2},$
(1)

where $n_{\rm b}$ and $\alpha$ are defined as $n_{\rm b}=n_{\rm n}+n_{\rm p}$
and $\alpha=(n_{\rm n}-n_{\rm p})/n_{\rm b}$. Here, $n_{\rm n}$ and $n_{\rm
p}$ are the neutron and the proton number density, respectively. The five
coefficients $w_{0}$, $n_{0}$, $S_{0}$, $K_{0}$ and $L$ are the nuclear
saturation parameters defining uniquely the EOS of the crust. The first three
parameters ($w_{0}$, $n_{0}$, and $S_{0}$) are well constrained from
terrestrial experiments (Oertel et al. 2017; Li et al. 2019), while the
remaining two parameters ($K_{0}$ and $L$) can be hardly constrained
experimentally. Still, there is progress in this direction and already certain
constraints can be set, predicting the following values: $K_{0}=240\pm 20$ MeV
(Shlomo et al. 2006) and $L=60\pm 20$ MeV (Li et al. 2019). In fact, the
magnetar QPOs in SGR 1806-20 and SGR 1900+14, if identified as crustal
torsional oscillations, provide some more stringent constraints to $L$, that
is $L_{\rm QPO}=58-73$ MeV (Sotani et al. 2018).

Table 1: The EOS parameters adopted in this study, where $\eta$ and $\varsigma$ are the combination of $K_{0}$ and $L$ given by $\eta=(K_{0}L^{2})^{1/3}$ and $\varsigma=(K_{0}^{4}L^{5})^{1/9}$. The corresponding transition density from spherical to cylindrical nuclei (SP–C) and that from cylindrical to slab-like nuclei (C–S) are also listed. The asterisk at the value of $K_{0}$ denotes the EOS model where cylindrical nuclei directly change to uniform matter.  $K_{0}$ (MeV) |  $L$ (MeV) |  $\eta$ (MeV) |  $\varsigma$ (MeV) | SP–C (fm-3) | C–S (fm-3) |
---|---|---|---|---|---|---
180 | 31.0 | $55.7$ | $67.8$ | 0.05887 | 0.07629 |
180 | 52.2 | $78.9$ | $90.5$ | 0.06000 | 0.07186 |
230 | 42.6 | $74.7$ | $90.1$ | 0.06238 | 0.07671 |
230 | 73.4 | $107$ | $122$ | 0.06421 | 0.07099 |
280 | 54.9 | $94.5$ | $113$ | 0.06638 | 0.07743 |
280∗ | 97.5 | $139$ | $156$ | 0.06678 | 0.06887 |
360 | 12.8 | $38.9$ | $56.4$ | 0.05777 | 0.08217 |
360 | 76.4 | $128$ | $152$ | 0.07239 | 0.07797 |

The OI-EOSs (Oyamatsu & Iida 2007) adopted in this study are constructed in
such a way that $w_{0}$, $n_{0}$, and $S_{0}$ are optimized for given values
of $K_{0}$ and $L$ to reproduce the experimental data for masses and charge
radii of stable nuclei, based on the extended Thomas-Fermi theory. The EOS
parameters adopted here are listed in Table 1. Therefore, for constructing the
crust equilibrium model, one needs to select two parameters for the EOS, i.e.,
$K_{0}$ and $L$, apart from the global parameters defining the star i.e., the
mass ($M$) and radius ($R$).

Torsional oscillations are characterized by the shear modulus, $\mu$, since
their frequency scales roughly as $f\propto\left(\mu/\rho\right)^{1/2}$ where
$\rho$ the local value of the density. For the phase composed of spherical
nuclei, the shear modulus is formulated as a function of the _ion number
density_ , the _charge number of nuclei_ , and the _Wigner-Seitz cell radius_
(Ogata & Ichimaru 1990; Strohmayer et al. 1991). When the phase is composed of
non-spherical (cylindrical and slab-like) nuclei, as proposed by Pethick &
Potekhin (1998), it behaves more like a liquid crystal. Then the shear modulus
in the phase of slab-like nuclei ceases to exist and at this part of the crust
torsional oscillations are not excited at least at a linear level. Thus,
torsional oscillations can be excited in two distinct areas of the crust a) in
the parts where the phase is composed of spherical and cylindrical nuclei, and
b) in the part where the phase is composed of cylindrical holes (tubes) and
spherical holes (bubbles) nuclei. This leads to double-layer torsional
oscillations in the crust, which can be visualized as a “lasagna sandwich”.
With such a structure, all QPOs observed in SGR 1806-20 and SGR 1900+14 can be
identified, see Sotani et al. (2019).

In this study, we focus on the newly observed QPOs in GRB 200415A. These QPOs
are of relatively higher frequency, and we suggest that they correspond to
overtones of crustal torsional oscillations. Thus, we systematically examine
the overtones of the torsional oscillations excited in a phase composed of
spherical and cylindrical nuclei. In addition, since some of the neutrons are
not bound in the nuclei inside the inner crust, one may take into account the
effect of such unbound neutrons on the torsional oscillations. In practice,
the ratio of the superfluid to the dripped neutron, $N_{\rm s}/N_{\rm d}$, in
the phase composed of spherical nuclei, is calculated by the band theory
(Chamel 2012). The corresponding ratio in the phase composed of cylindrical
nuclei is still unclear. So, in this study, we adopt the approach of Chamel
(2012) for spherical nuclei and $N_{\rm s}/N_{\rm d}$ as a parameter for
cylindrical nuclei, as in Sotani et al. (2018).

The governing equation for torsional oscillations is derived from the
linearized equations of motions (Schumaker & Thorne 1983), where we adopt the
relativistic Cowling approximation as in Sotani et al. (2007). Then, the
appropriate boundary conditions are imposed at the stellar surface and the
bottom of the crust (composed of cylindrical nuclei). In this way, a well-
defined eigenvalue problem is set for the frequency, $f$, of torsional
oscillations. For further details, we refer to Sotani et al. (2012, 2018).

Figure 1: The 1st overtone of the $\ell=2$ crustal torsional oscillations for
a neutron star model with $M=1.6M_{\odot}$, $R=12$km, and $N_{\rm s}/N_{\rm
d}=0$. The marks denote the torsional oscillation frequencies estimated for
the neutron star models constructed based on the crustal EOSs listed in Table
1, as a function of $\varsigma$. The thick-solid line stands for the quadratic
fitting given by Eq. (3).

The fundamental frequencies of crustal torsional oscillations scale roughly as
$f\sim v_{s}/R$, whereas the overtones scale as $f\sim v_{s}/\Delta R$, where
$v_{s}\sim\sqrt{\mu/\rho}$ is the shear velocity, and $\Delta R$ the crust
thickness (Hansen & Cioffi 1980). Since the crust thickness depends on both
the stellar compactness ($\Delta R/R\sim R/M$) and on the EOS parameters, as
shown in Sotani et al. (2017), the overtones generally depend on $K_{0}$ and
$L$. Actually, in Sotani et al. (2018) it was shown that the 1st overtone of
crustal torsional oscillations can be described as a function of the parameter

$\varsigma\equiv(K_{0}^{4}L^{5})^{1/9}.$ (2)

In fact, as it is shown in Fig. 1, the 1st overtone of a neutron star model
with $M=1.6M_{\odot}$, $R=12$ km, and $N_{\rm s}/N_{\rm d}=0$ is expressed as
a function of $\varsigma$. The colored marks correspond to the
eigenfrequencies, for some typical values of $K_{0}$, while the thick-solid
line represents the following quadratic fitting

${}_{\ell}t_{n}=d_{\ell n}^{(0)}+d_{\ell n}^{(1)}\varsigma_{100}+d_{\ell
n}^{(2)}\varsigma^{2}_{100}\,.$ (3)

Here, $\varsigma_{100}\equiv\varsigma/(100{\rm MeV})$, while the fitting
coefficients $d_{\ell n}$, depend on the values of the azimuthal quantum
number $\ell$, the nodal number in the corresponding eigenfunction $n$, the
stellar mass $M$, the radius $R$, and the ratio $N_{\rm s}/N_{\rm d}$.
Actually, in this study, we further confirm that not only the 1st overtone but
also the $n$-th overtones can be fitted by an equation of the form of Eq. (3).
Thus, hereafter, we will discuss the EOS dependence of the overtones using the
fitting formula given by Eq. (3) as a function of $\varsigma$. As discussed
earlier, the overtone frequencies generally increase with the stellar
compactness, $M/R$ (since $\Delta R/R\sim R/M$) and $f\sim v_{s}/\Delta R$ for
the overtones (Hansen & Cioffi 1980).

##  3 Comparison of the overtones with the observed QPOs

In this section, we make an attempt to identify which overtones of crustal
torsional oscillations match the QPO frequencies observed in GRB 200415A
keeping in mind that not all four QPOs were of the same significance. The
lower two were less prominent, which apart from underlying physical reasons
(which are unclear, due to the lack of knowledge of the excitation mechanism),
should be related also to the short duration of the signal ($\sim 3$ms).

As a first step, we will consider the correspondence between the 1st overtone
and the lowest QPO frequency extracted from GRB 200415A, i.e.,
835.9${}^{-84.7}_{+77.3}$ Hz (Castro-Tirado et al. 2021). In Fig. 2, we show
the 1st overtones for three characteristic neutron star models with
$(M,R)=(1.4M_{\odot},14{\rm km})$, $(1.6M_{\odot},12{\rm km})$, and
$(1.8M_{\odot},10{\rm km})$. The solid and dotted lines denote the frequencies
with $\ell=2$ and $\ell=10$, while the top and bottom panels correspond to the
results with $N_{\rm s}/N_{\rm d}=0$ and $N_{\rm s}/N_{\rm d}=1$,
respectively. From this figure, it is obvious that the dependencies of the
overtones on $\ell$ and $N_{\rm s}/N_{\rm d}$ is very weak, as mentioned in
Hansen & Cioffi (1980); Sotani et al. (2018). So, hereafter, we will only
consider the $\ell=2$ overtones for the neutron star model with $N_{\rm
s}/N_{\rm d}=0$. In addition, it is noticeable that the 1st overtones increase
with compactness as mentioned in the previous section. Actually, for the three
models considered here, $(M,R)=(1.4M_{\odot},14{\rm km})$,
$(1.6M_{\odot},12{\rm km})$ and $(1.8M_{\odot},10{\rm km})$, the compactness
is $M/R=0.148$, $0.197$ and $0.266$, correspondingly. As a result, for
identifying the 836 Hz QPO as the 1st overtone, the suitable values of
$\varsigma$ should increase as the neutron star compactness becomes higher.

Figure 2: The 1st overtones of the crustal torsional oscillations with
$\ell=2$ (solid line) and $\ell=10$ (dashed line) are shown as function of
$\varsigma$ given by Eq. (2) for three neutron star models with
$(M,R,M/R)=(1.4M_{\odot},14{\rm km},0.146)$, $(1.6M_{\odot},12{\rm km},0.197)$
and $(1.8M_{\odot},10{\rm km},0.266)$. The top panel corresponds to $N_{\rm
s}/N_{\rm d}=0$ while the bottom panel to $N_{\rm s}/N_{\rm d}=1$. In both
panels, for reference, we draw (with a shaded thick line) the lowest QPO
frequency found in GRB 200415A, i.e., 835.9${}^{-84.7}_{+77.3}$ Hz (Castro-
Tirado et al. 2021).  Figure 3: Suitable values of $\varsigma$ for identifying
the 835 Hz QPO as the 1st overtone for various stellar masses and radii with
$N_{\rm s}/N_{\rm d}=0$ are shown. The shaded region corresponds to the range
of $\varsigma=85.3-135.1$MeV, estimated for the fiducial values $L=60\pm 20$
and $K_{0}=240\pm 20$ MeV. The dashed lines define the range of
$\varsigma_{\rm QPO}=104.9-128.4$ MeV, corresponding to $L_{\rm QPO}=58-73$
MeV.

In Fig. 3, we show the suitable value of $\varsigma$ for identifying the 836
Hz QPO as the 1st overtone of crustal torsional oscillations for some typical
neutron star models. More specifically, the shaded region corresponds to
neutron star models for $\varsigma=85.3-135.1$ MeV, corresponding to
$K_{0}=240\pm 20$ MeV and $L=60\pm 20$ MeV. In addition, we draw the dashed
lines at $\varsigma=104.9$ and $\varsigma=128.4$ MeV, defining the estimated
range for the tighter range $L_{\rm QPO}=58-73$ MeV. The latter bounds have
been set once the magnetar QPOs in SGR 1806-20 and 1900+14 were identified as
crustal torsional oscillations (Sotani et al. 2018). In this figure, it is
noticeable that the suitable value of $\varsigma$, for identifying the 836 Hz
QPO as the 1st overtone, increases with $M$, for fixed neutron star radius,
and that some characteristic neutron star models, e.g., with masses
$M=1.4M_{\odot}$, $1.6M_{\odot}$, $1.8M_{\odot}$ and $R=10$ km, do not fit
with the observational data.

The next step is to associate all four QPO frequencies observed in GRB 200415A
with overtones of the crustal torsional oscillations. In fact, we find that it
is feasible to identify all four observed QPOs as the 1st, 2nd, 4th, and 10th
overtones for neutron stars models with certain values of mass and radius. In
Fig. 4, we show the suggested identification for one of the three neutron star
models considered earlier i.e., the one with $M=1.6M_{\odot}$ and $R=12$ km.
In this figure, the correspondence for $\varsigma=121.7$ MeV for the
aforementioned neutron star model is apparent. Since the overtones vary as
they depend on the neutron star mass and radius, the suitable value of
$\varsigma$ for identifying the observed QPO frequencies varies accordingly.

Figure 4: Correspondence of the observed QPOs to the overtones for the neutron
star model with $M=1.6M_{\odot}$ and $12$ km. The QPOs observed in GRB 200415A
can be identified with the 1st, 2nd, 4th, and 10th overtones of crustal
torsional oscillations, which leads to the suitable value of $\varsigma$ for
the adopted stellar model as $\varsigma=121.7$ MeV.

In Fig. 5, we show (by colored marks) the extracted values of $\varsigma$ for
a wider range of neutron star models. In all of them, we identified the
observed QPO frequencies as the 1st, 2nd, 4th, and 10th overtones. As in Fig.
3, in this figure we draw the fiducial range of $\varsigma$ (shaded region)
and the range for $\varsigma_{\rm QPO}$ (dashed lines) as constrained by QPO
observations. We note that for the stellar models, which are not considered in
Fig. 5, e.g., for a neutron star model with $M=1.8M_{\odot}$ and 10 km, the
QPO frequencies cannot be identified with the same set of overtones.

Figure 5: The suitable value of $\varsigma$ (colored marks) for simultaneous
identification of the four QPOs observed in GRB 200415A as the 1st, 2nd, 4th,
and 10th overtones of the crustal torsional oscillations for various neutron
star models. The meaning of the shaded region and dashed lines is the same as
in Fig. 3. The arrow denotes the particular model discussed in Fig. 4, where
the extracted value is $\varsigma=121.7$ MeV.

Fig. 6 is a mass vs radius diagram on which we extract the constraints set in
Fig. 5 for the mass, radius and $\varsigma$. The shaded region (the region
enclosed with the solid lines) corresponds to the mass vs radius constraints
obtained for $\varsigma=85.3-135.1$ MeV ($\varsigma_{\rm QPO}=104.9-128.4$
MeV).

Figure 6: Constraints on neutron star mass and radius set by the simultaneous
identification of all four observed QPOs in GRB 200415A as the 1st, 2nd, 4th,
and 10th overtones of the crustal torsional oscillations. Here we accepted
only the models meeting the constraints $\varsigma=85.3-135.1$ MeV (shaded
region) and $\varsigma_{\rm QPO}=104.9-128.4$ MeV (bounded by solid lines).

It has been shown, that one can estimate the mass and radius of low-mass
neutron stars, using another combination of $L$ and $K_{0}$, i.e.,
$\eta=(K_{0}L^{2})^{1/3}$ (Sotani et al. 2014, 2022)333A similar discussion is
also possible using another combination of the nuclear saturation parameter
given by $\eta_{\tau}=(-K_{\tau}L^{5})^{1/6}$, where $K_{\tau}$ is the isospin
dependence of incompressibility for asymmetric nuclear matter (Sotani & Ota
2022), or directly using the experimental observables, such as neutron skin
thickness or dipole polarizability for neutron-rich nuclei (Sotani & Naito
2023). According to this approach, for neutron star models with central
density, $\rho_{\rm c}\leq 2\rho_{0}$, its mass and gravitational redshift can
be expressed as a function of $\eta$ and the normalized $\rho_{\rm c}$. This
leads to a relation between the neutron star mass and radius, once $\eta$ is
fixed. Then we adopt a range for $\eta=70.6-118.5$ MeV, corresponding to
$\varsigma=85.3-135.1$ MeV and $\eta_{\rm QPO}=90.5-111.5$ MeV, corresponding
to $\varsigma_{\rm QPO}=104.9-128.4$ MeV. The expected region of the neutron
star mass and radius with $\rho_{\rm c}\leq 2\rho_{0}$, is given by the shaded
region and the region enclosed by the solid lines, on the bottom-right side in
Fig. 7. If the resultant region intersects with the region shown in Fig. 6,
one could set further constraints to the parameters of the neutron star model
corresponding to GRB 200415A. Unfortunately, the shaded area corresponding to
neutron star models with $\rho_{\rm c}\leq 2\rho_{0}$, does not intersect with
the region shown in Fig. 6. Thus, to narrow down the allowed region for the
neutron star models corresponding to GRB 200415A, we will make a further
assumption so that the expected mass and radius could somehow intersect with
the region shown in Fig. 6.

As seen in Fig. 7, the masses of the neutron star models with $\rho_{\rm
c}\leq 2\rho_{0}$ are quite small and they do not overlap with the constraint
mass vs radius area found in Fig. 6. Still, the neutron star radius hardly
changes at all in the domain of $M\simeq 0.5-1.5M_{\odot}$ for any EOSs. Based
on this observation we can extrapolate our mass vs radius domain as
constrained by the vertical dashed lines corresponding to $R=11.73$, $12.41$,
$13.03$, and $13.23$ km for $\eta=70.6$, $90.5$, $111.5$, and $118.5$ MeV. In
this way, we create an overlapping of the constrained region based on the
information extracted from GRB 200415A, as shown in Fig. 6, with the neutron
star models based on the allowed values of $\eta$ ($\eta_{\rm QPO}$).

Figure 7: The allowed range for the neutron star mass and radius obtained from
the QPOs observed in GRB 200415A, as shown in Fig. 6, is compared to the mass
and radius range derived via the empirical formula derived in Sotani et al.
(2014) (the right-bottom region). For defining the low mass models we used
$\eta=70.6-118.5$ MeV (shaded region) and $\eta_{\rm QPO}=90.5-111.5$ MeV (the
region enclosed by the solid lines). Here we made the assumption that
$\rho_{c}\leq 2\rho_{0}$. The (black and blue) dashed lines denote the neutron
star radii expected from the empirical relation for $\eta=70.6$, 90.5, 111.5,
and 118.5 MeV with $\rho_{c}=2\rho_{0}$.

Finally, the allowed mass vs radius area for the neutron star models
constrained by the observed QPOs in GRB 200415A is compared to the other
constraints on the neutron star mass and radius obtained from the astronomical
and experimental observations. Actually, in Fig. 8 we define, with the two
parallelograms, the allowed region for the neutron star model that matches the
information extracted from GRB 200415A. The outer (inner) parallelogram
corresponds to the expected region assuming $\eta=70.6-118.5$ MeV ($\eta_{\rm
QPO}=90.5-111.5$ MeV). Meanwhile, in the same figure, we plot the constraints
on the neutron star mass and radius set by (i) NICER observations for PSR
J0030+0451 and MSP J0740+6620 (Riley et al. 2019; Miller et al. 2019; Riley et
al. 2021; Miller et al. 2021); (ii) observations of x-ray bursts (Steiner et
al. 2013); and (iii) gravitational wave observations and constraints set by
GW170817. Especially, for GW170817, we draw the conservative constraint set by
the tidal deformability, predicting that the radius of a neutron star with
$M=1.4M_{\odot}$ should be $R_{1.4}\leq 13.6$km (Annala et al. 2018) and that
with $M=1.6M_{\odot}$ should be $R_{1.6}\geq 10.7$km (Bauswein et al. 2017),
together with more stringent constraints obtained from the combination of
multimessenger observations and nuclear theory, i.e.,
$R_{1.4}=11.0^{+0.9}_{-0.6}$ km (Capano et al. 2020) and
$R_{1.4}=11.75^{+0.86}_{-0.81}$ km (Dietrich et al. 2020).

Recently, a new upper neutron star mass limit was set by PSR J0952-0607, the
so-called “black widow”, that is $M=2.35\pm 0.17M_{\odot}$ (Romani et al.
2022) which is also shown in Fig. 8. In the high mass region we also drew the
limits set by causality (Lattimer 2012). Finally, the low-mass region enclosed
by the purple solid line is based on the constraints set for $L=60\pm 20$ and
$K_{0}=240\pm 20$ MeV, as in Fig. 7. For reference, in Fig. 8, we also show
the neutron star models constructed using realistic EOSs. Among the EOSs
adopted here, the Shen EOS has already been excluded from the gravitational
wave observation in GW170817, we keep it here as a reference since it is one
of the standard EOSs in astrophysics.

Figure 8: The neutron star mass and radius constrains in this study from GRB
200415A (double-parallelogram) are shown together with other observational and
experimental constraints (see the text for details). For reference, we also
plot the neutron star mass and radius theoretically constructed with several
realistic EOSs.

##  4 Conclusion

Quasi-periodic oscillations (QPOs) observed in the magnetar giant flares can
be an extremely useful tool for extracting neutron star properties. In this
study, we examine the possibility for identifying the QPOs of the recently
observed GRB 200415A as overtones of crustal torsional oscillations. As a
result, we find that the observed QPO frequencies can be identified as the
1st, 2nd, 4th, and 10th overtones with specific values of
$\varsigma=(K_{0}^{4}L^{5})^{1/9}$.

Then, by comparing the resultant value of $\varsigma$ for a variety of neutron
star models to the appropriate range of $\varsigma$ restricted by the values
of $L$ and $K_{0}$ determined by terrestrial experiments (or even the value of
$L$ constrained from the previous giant flare observations), we show the
possible range of values for the mass and radius of GRB 200415A. Furthermore,
assuming that the neutron star radius is almost the same as that of the
neutron star model, whose central density is twice the saturation density, we
derived a more stringent constraint for the neutron star mass and radius of
GRB 200415A. The outcome of our analysis is in good agreement, with the
constraints set by other observations in the gravitational and electromagnetic
spectrum (Fig. 8).

In this study, we worked under the assumption that the strength of the
magnetic field in GRB 200415A is
$\,\,\raisebox{-3.44444pt}{$\stackrel{{\scriptstyle\textstyle<}}{{\sim}}$
}10^{15}$G, so that torsional oscillations have a short damping time, but the
frequencies are still close to the values of pure crustal oscillations. If it
turns out that the magnetic field strength in GRB 200415A was higher than
$\sim 10^{15}$G, a revised analysis will be necessary, due to a significant
effect on the frequency spectrum of torsional oscillations. We emphasize that
the alternative model discussed in Castro-Tirado et al. (2021) is also a
viable solution.

## Acknowledgements

This work is supported in part by Japan Society for the Promotion of Science
(JSPS) KAKENHI Grant Numbers JP19KK0354 and JP21H01088, and by Pioneering
Program of RIKEN for Evolution of Matter in the Universe (r-EMU).

##  Appendix A Magnetic effects on the crustal torsional oscillations

In this study, we worked under the assumption that the strength of the
magnetic field in GRB 200415A is
$\,\,\raisebox{-3.44444pt}{$\stackrel{{\scriptstyle\textstyle<}}{{\sim}}$
}10^{15}$G, so that torsional oscillations have a short damping time, but the
frequencies are still close to the values of pure crustal oscillations. For
$B\geq 10^{15}$ G, there is a strong shift in the torsional oscillation
frequencies, and the damping due to the “continuous” spectrum is even
stronger, if one considers particular geometries for the magnetic field (Levin
2007; Colaiuda et al. 2009; van Hoven & Levin 2011; Colaiuda & Kokkotas 2011;
Gabler et al. 2012, 2018). For mixed poloidal-toroidal fields, the
magnetoelastic oscillations spectrum is becoming discrete (Colaiuda & Kokkotas
2012).

The effect of the magnetic field on pure torsional oscillations has been
studied in Sotani et al. (2007). There it was shown for a variety of neutron
star models and EOS that the shift in the torsional oscillations frequencies
due to magnetic field obeys the following formula

$\frac{{}_{\ell}f_{n}}{{}_{\ell}f^{(0)}_{n}}\approx\left[1+{}_{\ell}\alpha_{n}\left(\frac{B}{B_{\mu}}\right)^{2}\right]^{1/2},$
(4)

where ${}_{\ell}f^{(0)}_{n}$ is the torsional mode frequency of the $n$-th
overtone of a non-magnetized neutron star, while ${}_{\ell}f_{n}$ is the
frequency of the equivalent magnetized one with the same parameters ($M$, $R$,
EOS). $B$ is the strength of the surface magnetic field normalized by
$B_{\mu}=4\times 10^{15}\,{\rm G}$. This was the case also in the 2D linear
simulations in (Colaiuda & Kokkotas 2011, 2012) and it was almost confirmed in
the most recent 2D nonlinear studies with a poloidal magnetic field (Gabler et
al. 2018), where the oscillation frequencies are named as _magnetically
modified torsional modes_.

In the previous studies (Sotani et al. 2007), the coefficients
${}_{\ell}\alpha_{n}$, have been calculated only for $n=0,\,1$ and their
values vary from 0.3 to 0.5 for EOS NV (Negele & Vautherin 1973) and 0.4 to
1.5 for EOS DH (Douchin & Haensel 2001). Here we calculated ${}_{2}\alpha_{n}$
for larger values of $n$ as it is shown in Figure 9. The values seem to reach
a maximum value of about ${}_{2}\alpha_{n}\approx 2-2.5$ for EOS DH and
${}_{2}\alpha_{n}\approx 0.8-1.1$ for EOS NV. Thus the deviation of the
magnetized neutron star frequencies from those of the non-magnetized ones are
$\,\,\raisebox{-3.44444pt}{$\stackrel{{\scriptstyle\textstyle<}}{{\sim}}$
}3.4\%$ for the EOS NV and
$\,\,\raisebox{-3.44444pt}{$\stackrel{{\scriptstyle\textstyle<}}{{\sim}}$
}7.5\%$ for the EOS DH, if we assume $B\sim 10^{15}$ G. These values are still
within the limits of uncertainty ($\sim 10\%$) estimated in Castro-Tirado et
al. (2021). For values of the magnetic field significantly higher than $B\sim
10^{15}$ G, our approach would need to be modified.

Regarding the damping of crustal oscillations due to the Alfvén continuum in
the core, we remind that the 3.5ms duration of the high-frequency QPOs
observed in GRB 200415A is consistent with the example of a timescale of $\sim
5$ms for the damping of a crustal $n=1$ shear mode in a model with a magnetic
field strength of $\sim 10^{15}$G shown in Gabler et al. (2018), which
justifies our assumptions.

Figure 9: The coefficients of ${}_{2}a_{n}$ in Eq. (4) for the neutron star
models constructed with two different crust EOSs, NV (Negele & Vautherin 1973)
and DH (Douchin & Haensel 2001). For the interior (core) we have used the APR
and WFF3 EOS. The open and filled marks correspond to results with NV and DH
EOS, while the solid, dashed, and dotted lines correspond to the results for
the $1.4M_{\odot}$, $2.0M_{\odot}$ neutron star models for the APR EOS, and
$1.4M_{\odot}$ neutron star models for the WFF3 EOS.

##  Appendix B Can the 836 Hz QPO be the 2nd overtone?

In the main part of the article, we identified the lowest QPO frequency at 836
Hz as the 1st overtone. In this appendix, we examine the possibility of
identifying the lowest QPO as the 2nd overtone. In Fig. 10 we plot the 1st and
2nd overtone of crustal torsional oscillations for a neutron star model with
$1.4M_{\odot}$ and $14$ km as a function of $\varsigma$. In the same figure we
draw the observed QPO frequency and the aforementioned values for $\varsigma$
and $\varsigma_{\rm QPO}$. It is apparent that the 836 Hz QPO frequency can be
identified as the 2nd overtone. In fact, as it is shown in Fig. 11, all four
QPOs observed in GRB 200415A can be identified as the 2nd, 5th, 8th, and 16th
overtones, if $\varsigma=142.1$ MeV. However, as it was mentioned earlier, the
overtone frequencies increase with compactness, which implies that the
suitable $\varsigma$ values for identifying the observed QPOs should be
larger, beyond the accepted values for $\varsigma$. Therefore, for being in
agreement with the range $\varsigma=85.3-135.1$ MeV, we would need to consider
neutron star models with very small and somehow unphysical compactness.

Figure 10: The 1st (dotted line) and 2nd (solid line) overtones of the
$\ell=2$ crustal torsional oscillations for a neutron star model with
$(M,R)=(1.4M_{\odot},14{\rm km})$ and $N_{\rm s}/N_{\rm d}=0$ is shown as
functions of $\varsigma$. For reference, we plot the lowest observed QPO
frequency (red shaded area) in GRB 200415A, i.e., 835.9${}^{-84.7}_{+77.3}$ Hz
(Castro-Tirado et al. 2021). The ranges for $\varsigma$ (dark shaded area) and
$\varsigma_{\rm QPO}$ (the area between the dashed lines) are also drawn for
comparison.  Figure 11: Identification of the QPOs observed in GRB 200415A as
crustal torsional oscillation overtones for a neutron star model with
$M=1.4M_{\odot}$ and $R=14$ km, assuming that the lowest QPO in GRB 200415A
corresponds to the 2nd overtone.

## References

  * Abbott et al. (2017) Abbott, B. P. et al. (LIGO Scientific Collaboration and Virgo Collaboration), Phys. Rev. Lett., 119, 161101
  * Abbott et al. (2019) Abbott, B. P. et al. (LIGO Scientific Collaboration and Virgo Collaboration), Phys. Rev. X, 9, 011001
  * Andersson & Kokkotas (1996) Andersson N., & Kokkotas K. D. 1996, Phys. Rev. Lett., 677, 4134
  * Andersson & Kokkotas (1998) Andersson N., & Kokkotas K. D. 1998, MNRAS, 299, 1059
  * Annala et al. (2018) Annala E., Gorda T., Kurkela A., & Vuorinen A. 2018, Phys. Rev. Lett., 120, 172703
  * Antoniadis et al. (2013) Antoniadis J., Freire P. C. C., Wex N., et al. 2013, Science, 340, 6131
  * Bauswein et al. (2017) A. Bauswein, O. Just, H.-T. Janka, and N. Stergioulas, 2017, ApJL, 850, L34
  * Capano et al. (2020) Capano C. D., Tews I., Brown S. M., et al. 2020, Nature Astron. 4, 625
  * Castro-Tirado et al. (2021) Castro-Tirado A. J., Østgaard N., Göğüş E., et al. 2021, Nature, 600, 621
  * Chamel (2012) Chamel N. 2012, Phys. Rev. C, 85, 035801
  * Chatziioannou (2020) Chatziioannou, K., Gen. Rel. Grav., 52, 109
  * Dietrich, Hinderer & Samajdar (2021) Dietrich, T., Hinderer, T. and Samajdar, A., Gen. Rel. Grav., 53, 27
  * Colaiuda et al. (2009) Colaiuda A., Beyer, H., & Kokkotas K. D. 2009, MNRAS, 396, 1441
  * Colaiuda & Kokkotas (2011) Colaiuda A., & Kokkotas K. D. 2011, MNRAS, 414, 3014
  * Colaiuda & Kokkotas (2012) Colaiuda A., & Kokkotas K. D. 2012, MNRAS, 423, 811
  * Cromartie et al. (2020) Cromartie H. T., Fonseca E., Ransom S. M., et al. 2020, Nature Astronomy, 4, 72
  * Demorest et al. (2010) Demorest P., Pennucci T., Ransom S., et al. 2010, Nature, 467, 1081
  * Dietrich et al. (2020) Dietrich T., Coughlin M. W., Pang P. T. H., et al. 2020, Science 370, 1450
  * Doneva et al. (2013) Doneva D. D., Gaertig E., Kokkotas K. D., & Krüger C. 2013, Phys. Rev. D, 88, 044052
  * Douchin & Haensel (2001) Douchin F., & Haensel P. 2001, Astron. Astrophys., 380, 151
  * Gabler et al. (2011) Gabler M., Cerdá Durán P., Font J. A., et al. 2011, MNRAS, 410, L37
  * Gabler et al. (2012) Gabler M., Cerdá Durán P., Stergioulas N., et al. 2012, MNRAS, 421, 2054
  * Gabler et al. (2013) Gabler M., Cerdá Durán P., Stergioulas N., et al. 2013, Phys. Rev. Let., 111, 211102
  * Gabler et al. (2016) Gabler M., Cerdá Durán P., Stergioulas N., et al. 2016, MNRAS, 460, 4242-4257
  * Gabler et al. (2018) Gabler M., Cerdá Durán P., Stergioulas N., et al. 2018, MNRAS, 476, 4199
  * Gearheart et al. (2011) Gearheart M., Newton W. G., Hooker J., & Li B. A. 2011, MNRAS, 418, 2343
  * Haensel et al. (2006) Haensel P., Potekhin A. Y., & Yakovlev D. G., in Neutron Stars 1: Equation of State and Structure (Springer, 2006)
  * Hansen & Cioffi (1980) Hansen C., & Cioffi D. F. 1980, ApJ, 238, 740
  * Israel et al. (2005) Israel G. et al., 2005 Astrophys. J., 628, L53
  * Krüger et al. (2021) Krüger C. J., Kokkotas K. D., Manoharan P., & Völkel S. H. 2021, Front. Astron. Space Sci. 8, 736918
  * Lattimer (2012) Lattimer J. M. 2012, Ann. Rev. Nucl. Part. Sci. 62, 485
  * Leahy & Li (1995) Leahy D. A., & Li L. 1995, MNRAS, 277, 1177
  * Levin (2007) Levin Y. 2007, MNRAS, 377, 159
  * Li et al. (2019) Li B. A., Krastev P. G., Wen D. H., & Zhang N. B. 2019, Euro. Phys. J. A, 55, 117
  * Mahlmann et al. (2019) Mahlmann, J. F., Akgün, T., Pons, J. A., Aloy, M. A. & Cerdá-Durán, P. 2019, MNRAS, 490, 4858
  * Miller et al. (2019) Miller M. C., Lamb F. K., Dittmann A. J., et al. 2019, ApJL, 887, L24
  * Miller et al. (2021) Miller M. C., Lamb F. K., Dittmann A. J., et al., 2021, ApJL, 918, L28
  * Morozova et al. (2018) Morozova V., Radice D., Burrows A., & Vartanyan D. 2018, ApJ, 861, 10
  * Negele & Vautherin (1973) Negele J. W., & Vautherin D. 1973, Nucl. Phys., A207, 298 (NV)
  * Oertel et al. (2017) Oertel M., Hempel M., Klähn T., & Typel S. 2017, Rev. Mod. Phys., 89, 015007.
  * Ogata & Ichimaru (1990) Ogata S., & Ichimaru S. 1990, Phys. Rev. A, 42, 4867
  * Oyamatsu & Iida (2003) Oyamatsu K., & Iida K., 2003 Prog. Theor. Phys., 109, 631
  * Oyamatsu & Iida (2007) Oyamatsu K., & Iida K. 2007, Phys. Rev. C, 75, 015801
  * Passamonti & Andersson (2012) Passamonti A., & Andersson N. 2012, MNRAS, 419, 638
  * Passamonti & Lander (2014) Passamonti A., Lander S. K., MNRAS, 438, 156
  * Pechenick et al. (1983) Pechenick K. R., Ftaclas C., & Cohen J. M. 1983, ApJ, 274, 846.
  * Pethick & Potekhin (1998) Pethick C. J., & Potekhin A. Y. 1998, Phys. Lett. B427, 7
  * Poutanen & Gierlinski (2003) Poutanen J., & Gierlinski M. 2003, MNRAS, 343, 1301
  * Psaltis & Özel (2014) Psaltis D., & Özel F. 2014, ApJ, 792, 87
  * Riley et al. (2019) Riley T. E., Watts A. L., Bogdanov S., et al. 2019, ApJL, 887, L21
  * Riley et al. (2021) Riley T. E., Watts A. L., Ray P. S., et al. 2021, ApJL, 918, L27
  * Romani et al. (2022) Romani R. W., Kandel D., Filippenko A. V., et al. 2022, ApJL, 934, L17
  * Schumaker & Thorne (1983) Schumaker B. L., & Thorne K. S. 1983, MNRAS, 203, 457
  * Shlomo et al. (2006) Shlomo S., Kolomietz V. M., & Colò G. 2006, Eur. Phys. J. A, 30, 23
  * Sotani & Miyamoto (2018) Sotani H., & Miyamoto U. 2018, Phys. Rev. D, 98, 044017
  * Sotani & Ota (2022) Sotani H., & Ota S. 2022, Phys. Rev. D, 106, 103005
  * Sotani & Naito (2023) Sotani H., & Naito T. 2023, Phys. Rev. C, in press
  * Sotani et al. (2004) Sotani H., Kohri K., & Harada T. 2004, Phys. Rev. D, 69, 084008
  * Sotani et al. (2007) Sotani H., Kokkotas K. D., & Stergioulas N. 2007, MNRAS, 375, 261
  * Sotani et al. (2011) Sotani H., Yasutake N., Maruyama T., & Tatsumi T. 2011, Phys. Rev. D, 83, 024014
  * Sotani et al. (2012) Sotani H., Nakazato K., Iida K., & Oyamatsu K. 2012, Phys. Rev. Lett., 108, 201101
  * Sotani et al. (2013) Sotani H., Nakazato K., Iida K., & Oyamatsu K. 2013, MNRAS, 428, L21
  * Sotani et al. (2014) Sotani H., Iida K., Oyamatsu K., & Ohnishi A. 2014, Prog. Theor. Exp. Phys. 2014, 051E01
  * Sotani et al. (2017) Sotani H., Iida K., & Oyamatsu K. 2017, MNRAS, 470, 4397
  * Sotani et al. (2018) Sotani H., Iida K., & Oyamatsu K. 2018, MNRAS, 479, 4735
  * Sotani et al. (2019) Sotani H., Iida K., & Oyamatsu K. 2019, MNRAS, 489, 3022
  * Sotani et al. (2021) Sotani H., Takiwaki T., & Togashi H. 2021, Phys. Rev. D, 104, 123009
  * Sotani et al. (2022) Sotani H., Nishimura N., & Naito T. 2022, Prog. Theor. Exp. Phys. 2022, 041D01
  * Steiner & Watts (2009) Steiner A. W., & Watts A. L. 2009, Phys. Rev. Lett., 103, 181101
  * Steiner et al. (2013) Steiner A. W., Lattimer J. M., & Brown E. F. 2013, ApJL, 765, L5
  * Strohmayer & Watts (2005) Strohmayer T. E., & Watts A. L. 2005, ApJ, 632, L111
  * Strohmayer & Watts (2006) Strohmayer T. E., & Watts A. L. 2006, ApJ, 653, 593
  * Strohmayer et al. (1991) Strohmayer T., van Horn H. M., Ogata S., et al. 1991, ApJ, 375, 679
  * Sun et al. (2023) Sun, X., Miao, Z., Sun, B. & Li, A., ApJ, 942, 55
  * Suvorov et al. (2022) Suvorov A. G., Kuan H. J., & Kokkotas K. D. 2022 A&A, 664, A177
  * van Hoven & Levin (2011) van Hoven M., & Levin Y. 2011, MNRAS, 410, 1036
  * Torres-Forné et al. (2019) Torres-Forné A., Cerdá-Durán P., Passamonti A., et al. 2019, MNRAS, 482, 3967
  * Xiao et al. (2022) Xiao S. et al., 2022, arXiv:2205.02186.