# Determining the equation of state of neutron stars with Einstein Telescope
using tidal effects and r-mode excitations from a population of binary
inspirals

Pawan Kumar Gupta1,2 Anna Puecher1,2 Peter T.H. Pang1,2 Justin Janquart1,2
Gideon Koekoek1,3 Chris Van Den Broeck1,2 1Nikhef – National Institute for
Subatomic Physics, Science Park 105, 1098 XG Amsterdam, The Netherlands
2Institute for Gravitational and Subatomic Physics (GRASP), Utrecht
University, Princetonplein 1, 3584 CC Utrecht, The Netherlands  3Department of
Gravitational Waves and Fundamental Physics, Maastricht University, P.O. Box
616, 6200 MD Maastricht, The Netherlands

###### Abstract

Third-generation gravitational wave (GW) observatories such as Einstein
Telescope (ET) and Cosmic Explorer (CE) will be ideal instruments to probe the
structure of neutron stars through the GWs they emit when undergoing binary
coalescence. In this work we make predictions about how well ET in particular
will enable us to reconstruct the neutron star equation of state through
observations of tens of binary neutron star coalescences with signal-to-noise
ratios in the hundreds. We restrict ourselves to information that can be
extracted from the inspiral, which includes tidal effects and possibly r-mode
resonances. In treating the latter we go beyond the Newtonian approximation,
introducing and utilizing new universal relations. We find that the ability to
observe resonant r-modes would have a noticeable impact on neutron star
equation of state measurements with ET.

##  I Introduction

The detection with Advanced LIGO Aasi  _et al._ (2015) and Advanced Virgo
Acernese  _et al._ (2015) of gravitational waves (GWs) from the binary neutron
star (BNS) coalescences GW170817 Abbott  _et al._ (2017a) and GW190425 Abbott
_et al._ (2020a), together with electromagnetic observations Soares-Santos
_et al._ (2017); Cowperthwaite  _et al._ (2017); Abbott  _et al._ (2017b) has
already had a significant impact on our insight into the structure of neutron
stars; for a recent review, see e.g. Dietrich  _et al._ (2021). Even so, the
neutron star equation of state (EOS) remains poorly constrained. This is
expected to change with the advent of third-generation GW detectors such as
Einstein Telescope (ET) Punturo  _et al._ (2010); Hild  _et al._ (2011) and
Cosmic Explorer (CE) Abbott  _et al._ (2017c); Reitze  _et al._ (2019), which
are likely to see hundreds of thousands of binary neutron star coalescences,
of which hundreds may have signal-to-noise ratios (SNRs) in excess of 100; see
Samajdar  _et al._ (2021) for recent estimates. Whereas with current detectors
we only have access to the inspiral signal, ET and CE will also probe the
post-merger Maggiore  _et al._ (2020); Evans  _et al._ (2021).

In this paper we will quantify to what extent the equation of state –
essentially pressure versus density, $p(\rho)$ – of neutron stars can be
probed with third-generation GW observatories, based on inspiral physics
alone. The main EOS-related effect that enters the inspiral signal is that of
tides on the neutron stars Flanagan and Hinderer (2008); Hinderer (2008).
These can already be investigated with second-generation detectors Del Pozzo
_et al._ (2013); Agathos  _et al._ (2015); Lackey and Wade (2015), as has
indeed been done with GW170817 Abbott  _et al._ (2017a, 2018, 2020b). Tidal
effects are the most noticeable at high frequencies. However, given
sufficiently high SNR, it may also be possible to see resonant r-modes
Flanagan and Racine (2007), gravitomagnetic excitations of Rossby modes that
are induced when the monotonically increasing GW frequency reaches an
associated resonance frequency of one of the neutron stars. Assuming slowly
rotating neutron stars, these will mainly manifest themselves at lower
frequencies, i.e. tens of Hz. Though the imprints of resonant r-modes on GW
signals will generally be weak, they can be within reach of third-generation
observatories Poisson (2020). Moreover, as recently indicated by Ma et al. Ma
_et al._ (2020), EOS measurements may gain from the ability to access r-mode
information in addition to tidal effects. In this study we will focus on ET,
in part because it is predicted to have excellent low frequency sensitivity.

In Ma  _et al._ (2020), Fisher matrix estimates were made to assess the
measurability of tidal deformabilities and other parameters related to neutron
stars in binary inspiral, with or without resonant r-modes. For the purpose of
initial estimates, in the latter paper the r-modes were largely treated in the
Newtonian approximation, but it is known that relativistic corrections on both
the resonance frequency and the size of the induced GW phase shift can be
sizeable; see e.g. Idrisy  _et al._ (2015); Gupta  _et al._ (2021); Jiménez
Forteza  _et al._ (2018). Here we do take such corrections into account, and
we introduce and use new universal relations for both. Moreover, we want to
know what is the impact of having a sizeable number of high-SNR BNS signals at
one’s disposal. Ideally we would want to perform Bayesian parameter estimation
simulations, but limited computational resources prevent us from doing so. As
a compromise, we will still use the Fisher matrix to obtain multivariate
Gaussian approximations for the likelihoods of individual sources, but we
sample from those to obtain estimates for EOS measurements with input from
multiple detected sources. Recently, catalogs of binary black hole and binary
neutron star detections by ET and CE were simulated Samajdar  _et al._ (2021)
based on predictions for merger rate as a function of redshift; we will draw
from this work to obtain a representative sample for the 20 loudest BNS
signals that ET is likely to detect based on its predicted sensitivity curve.
Based on these we will arrive at estimates for how accurately one will be able
to reconstruct $p(\rho)$ with observations by ET, with or without r-modes.

The rest of this paper is structured as follows. In Sec. II we describe our
treatment of resonant r-modes in terms of universal relations. In Sec. III we
explain the setup of our analyses. Results are given in Sec. IV. Finally, an
overview and conclusions are presented in Sec. V. Unless stated otherwise,
throughout this paper we set $G=c=1$.

##  II Resonant r-modes: Frequencies, phase shifts, and universal relations

As two neutron stars spiral towards each other, the GW frequency increases
monotonically, and at one or more points in time it can become equal to an
internal resonance frequency of one of the neutron stars. The resulting
excitation takes away part of the orbital energy, which speeds up the orbital
motion; this in turn gets imprinted upon the phasing of the GW signal.
Assuming the two neutron star each undergo near-instantaneous resonances at
respective frequencies $f_{0}^{(1,2)}$, the resulting change in the frequency-
domain GW phase with respect to the point particle case can be modeled as
Flanagan and Racine (2007); Flanagan and Hinderer (2008)

$\displaystyle\Psi_{r}$ $\displaystyle=$
$\displaystyle\left(1-\frac{f}{f_{0}^{(1)}}\right)\Delta\Phi_{1}\Theta(f-f_{0}^{(1)})$
(1)
$\displaystyle+\left(1-\frac{f}{f_{0}^{(2)}}\right)\Delta\Phi_{2}\Theta(f-f_{0}^{(2)}),$

where $\Theta(x)$ denotes the usual step function. Specializing to r-mode
excitations, the shifts $\Delta\Phi_{i}$, $i=1,2$ take the form Ma  _et al._
(2020)

$\Delta\Phi_{i}=-2\frac{5\pi^{2}}{192}\left(\frac{4}{3}\right)^{2/3}\frac{\Omega^{2/3}_{i}}{\mathcal{M}^{10/3}}\mathcal{I}_{i}.$
(2)

Here $\Omega_{i}$, $i=1,2$ are the angular rotation frequencies of the neutron
stars, and $\mathcal{M}=(m_{1}m_{2})^{3/5}/(m_{1}+m_{2})^{1/5}$ is the chirp
mass associated with the component masses $m_{1}$, $m_{2}$. The r-mode
couplings $\mathcal{I}_{i}$ take the form

$\mathcal{I}_{i}=(\bar{I}^{r}_{i})^{2}\,m_{i}^{4}\,\text{sin}^{2}(\psi_{i})\,\text{cos}^{4}(\psi_{i}/2)\,\left(1-\frac{m_{i}}{M}\right),$
(3)

where $M=m_{1}+m_{2}$ is the total mass, and $\psi_{i}$ are the angles between
the neutron stars’ spin vectors $\mathbf{S}_{i}$ and the orbital angular
momentum. The r-mode _overlaps_ $\bar{I}^{r}_{i}$ depend on the neutron stars’
internal structure.

Let us now discuss in turn our treatment of r-mode frequencies and overlaps.

###  II.1 Universal relation for the r-mode frequencies

For a slowly rotating neutron stars and in the Newtonian limit, the r-mode
frequencies in the co-rotating frame are proportional to the angular rotation
frequency $\Omega$, and given by Flanagan and Racine (2007)

$\omega_{\ell m}=-\frac{2m}{\ell(\ell+1)}\Omega.$ (4)

For given values of the spherical harmonic indices $(\ell,m)$, the
corresponding resonance frequency in the inertial frame is

$\omega_{0}=m\Omega+\omega_{\ell m}=m\Omega-\frac{2m}{\ell(\ell+1)}\Omega.$
(5)

The associated GW frequencies appearing in Eq. (1) are then
$f_{0}=2\omega_{0}/(2\pi)$, or

$f_{0}=\frac{1}{\pi}(m-\kappa)\,\Omega,$ (6)

with $\kappa=2m/\ell(\ell+1)$ in the Newtonian limit. For simplicity, in this
work we will assume barotropic neutron stars, so that $|m|=l$, and focus on
$(\ell,m)=(2,2)$; within the Newtonian framework one then has $\kappa=2/3$.

As shown in Idrisy  _et al._ (2015), in a relativistic treatment of slowly
rotating stars, the resonance frequency can differ significantly from the
Newtonian value. In that work, numerical values for $\kappa$ were computed for
relativistic neutron stars described by a variety of tabulated EOSs, and a
fitting formula was obtained for the function $\kappa(C)$, with $C=m/R$ the
compactness, where $R$ is a neutron star’s radius. However, in this work we
are interested in the imprint of resonant r-modes on GW emission, where it is
more convenient to have $\kappa$ as function of the dimensionless neutron star
tidal deformability $\Lambda$, since this is a parameter that enters directly
into the waveform.

With this in mind, for each value of $\kappa$ reported in Idrisy  _et al._
(2015) (Table II), we computed $\Lambda$ for each EOS and compactness listed.
This was done using the TOV solver available in LALSuite LIGO Scientific
Collaboration (2018). As seen in Fig. 1, there is a clear dependence of
$\kappa$ on the logarithm of $\Lambda$. Using a least-squares fitting method,
we find that the functional dependence is approximated well by the universal
relation

$\kappa=0.3668+0.0498\log(\Lambda)-0.0025\log^{2}(\Lambda),$ (7)

with fitting residuals at the $\mathcal{O}(1\%)$ level.

Figure 1: Top panel: values of $\kappa$ and corresponding tidal
deformabilities $\Lambda$ computed for a variety of EOSs and compactnesses
listed in Idrisy  _et al._ (2015) (green dots), and the fit for
$\kappa(\Lambda)$ of Eq. 7. Bottom panel: fitting residuals.

###  II.2 Universal relation for the r-mode overlap

Next we turn to the r-mode overlap. In the Newtonian approximation one has Ma
_et al._ (2020)

$\bar{I}^{r}=\sqrt{\frac{1}{m^{5}}\int_{0}^{R}\rho r^{6}dr},$ (8)

where $\rho$ is the density of the neutron star, $m$ its mass, and $R$ its
radius. In Ma  _et al._ (2020), a universal relation was obtained for the
Newtonian $\bar{I}^{r}$ as function of $\Lambda$. In this work we instead
follow the relativistic description of Gupta  _et al._ (2021), where the
following expression was found:

$\bar{I}^{r}=\sqrt{\frac{15}{4\pi}(\Sigma_{\text{stat}}-\Sigma_{\text{irr}})}.$
(9)

Here $\Sigma_{\text{stat}}$ and $\Sigma_{\text{irr}}$ are dimensionless static
and irrotational magnetic Love numbers, respectively. Universal relations for
these quantities as function of $\Lambda$ were derived in Jiménez Forteza  _et
al._ (2018), which for completeness we display again here:

$\displaystyle\log(-\Sigma_{\text{irr}})$
$\displaystyle=-2.03+0.487\log(\Lambda)+9.69\times 10^{-3}\log^{2}(\Lambda)$
$\displaystyle\,\,\,\,\,\,+1.03\times 10^{-3}\log^{3}(\Lambda)-9.37\times
10^{-5}\log^{4}(\Lambda)$ $\displaystyle\,\,\,\,\,\,+2.24\times
10^{-6}\log^{5}(\Lambda),$ (10) $\displaystyle\log(\Sigma_{\text{stat}})$
$\displaystyle=-2.66+0.786\log(\Lambda)-1.00\times 10^{-2}\log^{2}(\Lambda)$
$\displaystyle\,\,\,\,\,\,+1.28\times 10^{-3}\log^{3}(\Lambda)-6.37\times
10^{-5}\log^{4}(\Lambda)$ $\displaystyle\,\,\,\,\,\,+1.18\times
10^{-6}\log^{5}(\Lambda).$ (11)

Together with Eq. (9), these yield a relation for $\bar{I}^{r}(\Lambda)$.

##  III Setup of the analyses

In this section we describe the waveform model used, the Fisher matrix
formalism for obtaining a multivariate Gaussian approximation of likelihoods,
the parameterization used for reconstructing EOSs, and our framework for
performing such reconstructions using multiple GW detections.

###  III.1 Waveform model

We will focus on the inspiral, and use an analytic, frequency domain waveform
following the stationary phase approximation Sathyaprakash and Dhurandhar
(1991), which takes the general form

$\displaystyle\tilde{h}(f)=\mathcal{A}f^{-7/6}e^{i\Psi(f)}.$ (12)

Here $\mathcal{A}$ collects parameters appearing in the amplitude: masses, sky
position, orientation of the orbital plane, and luminosity distance to the
source. The phase $\Psi(f)$ can be written as

$\Psi(f)=\Psi_{\rm PP}(f)+\Psi_{\rm SO}(f)+\Psi_{\rm T}(f)+\Psi_{r}(f),$ (13)

where $\Psi_{\rm PP}$ has point particle contributions up to 3.5 post-
Newtonian (PN) order, and $\Psi_{\rm SO}$ contains spin-orbit effects at
1.5PN; we consider spin-spin effects to be negligible in the case of binary
neutron stars. (For a review of the post-Newtonian approximation, see Blanchet
(2002).) $\Psi_{\rm T}$ contains tidal effects at 5PN and 6PN orders Vines
_et al._ (2011). Finally, $\Psi_{r}$ is as given in Eq. (1), and we note that
to good approximation, the angles $\psi_{1,2}$ entering Eq. (3) can usually be
considered approximately constant even in the presence of spin precession
Kidder  _et al._ (1993); Kidder (1995); Apostolatos  _et al._ (1994).

Without the r-mode contribution, the phase depends on the 8 free parameters

$\displaystyle\vec{\theta}_{nr}=({t_{c},\phi_{c},m_{1},m_{2},\tilde{\Lambda},\delta\tilde{\Lambda},\chi_{sz},\chi_{az}}),$
(14)

whereas with r-modes included, the free parameters are

$\displaystyle\vec{\theta}_{r}=({t_{c},\phi_{c},m_{1},m_{2},\tilde{\Lambda},\delta\tilde{\Lambda},\Omega_{1},\Omega_{2},\psi_{1},\psi_{2}}).$
(15)

Here $t_{c}$ and $\phi_{c}$ are respectively the time and phase of
coalescence; $m_{i}$, $i=1,2$ the component masses;
$\chi_{sz}=(\chi_{1z}+\chi_{2z})/2$ and $\chi_{az}=(\chi_{1z}-\chi_{2z})/2$
respectively the symmetric and antisymmetric dimensionless spins, with
$\chi_{1z}$ and $\chi_{2z}$ the spin components in the direction of orbital
angular momentum; $\Omega_{i}$ the spin angular frequencies of the neutron
stars; and $\psi_{i}$ the angles between the dimensionless spins and the
orbital angular momentum. We note that in the detector frame, the component
masses $m_{1,2}^{\rm det}$ that directly enter the waveform are redshifted
with respect to the source frame masses $m_{1,2}$: one has $m_{1,2}^{\rm
det}=(1+z)\,m_{1,2}$, with $z$ the redshift. With a global network of third-
generation (3G) observatories, for BNS signals the luminosity distance can be
measured to $\mathcal{O}(1\%)$ accuracy for the loudest sources Zhao and Wen
(2018), and with a cosmological model (e.g. from Planck Ade  _et al._ (2014)),
distance can be converted to redshift. We will assume that other 3G
observatories will be operational at the same time as ET, so that for the
high-SNR sources considered below, the distance uncertainties will be
sufficiently small that uncertainties on redshift can be neglected in the
conversion between source frame and detector frame masses. Finally,
$\tilde{\Lambda}$ and $\delta\tilde{\Lambda}$ are related to the individual
neutron stars’ tidal deformability parameters $(\Lambda_{1},\Lambda_{2})$
through Wade  _et al._ (2014)

$\displaystyle\tilde{\Lambda}$ $\displaystyle=$
$\displaystyle\frac{8}{3}\left[(1+7\eta-31\eta^{2})(\Lambda_{1}+\Lambda_{2})\right.$
(16)
$\displaystyle\,\,\,\,\,\,\,+\left.\sqrt{1-4\eta}(1+9\eta-11\eta^{2})(\Lambda_{1}-\Lambda_{2})\right],$
$\displaystyle\delta\tilde{\Lambda}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\left[\sqrt{1-4\eta}(1-\frac{13272}{1319}\eta+\frac{8944}{1319}\eta^{2})(\Lambda_{1}+\Lambda_{2})\right.$
$\displaystyle\,\,\,\,\,\,\,+\left.(1-\frac{15910}{1319}\eta+\frac{32850}{1319}\eta^{2}+\frac{3380}{1319}\eta^{3})(\Lambda_{1}-\Lambda_{2})\right],$

where $\eta=m_{1}m_{2}/(m_{1}+m_{2})^{2}$ is the symmetric mass ratio.

When comparing the cases with and without r-modes, we can use that

$\chi_{iz}=\bar{I}_{i}\Omega_{i}m_{i}\text{cos}(\psi_{i})$ (18)

for $i=1,2$. Here $\bar{I}_{i}$ are normalized moments of inertia, for which
one has the universal relation Yagi and Yunes (2013)

$\displaystyle\log\bar{I}$ $\displaystyle=$ $\displaystyle 1.47+8.17\times
10^{-2}\log(\Lambda)+1.49\times 10^{-2}\log^{2}(\Lambda)$
$\displaystyle+2.87\times 10^{-4}\log^{3}(\Lambda)-3.64\times
10^{-5}\log^{4}(\Lambda).$

###  III.2 Spectral parametrization of the EOS

For the EOS we will use the so-called spectral parameterization in terms of
the adiabatic index $\Gamma(p)$, defined as Koliogiannis and Moustakidis
(2019); Lindblom (2010)

$\Gamma(p)=\frac{\epsilon+p}{p}\frac{dp}{d\epsilon},$ (20)

where $\epsilon$ is energy density and $p$ is pressure. The EOS $\epsilon(p)$
is obtained from the adiabatic index by writing the above equation as

$\frac{d\epsilon}{dp}=\frac{\epsilon(p)+p}{p\Gamma(p)},$ (21)

or

$\epsilon(p)=\frac{\epsilon_{0}}{\mu(p)}+\frac{1}{\mu(p)}\int_{p_{0}}^{p}\frac{\mu(p^{\prime})}{\Gamma(p^{\prime})}dp^{\prime},$
(22)

where

$\mu(p)=\exp\left(-\int_{p_{0}}^{p}\frac{1}{p^{\prime}\Gamma(p^{\prime})}dp^{\prime}\right),$
(23)

with $\epsilon_{0}=\epsilon(p_{0})$ a constant of integration. The EOS
$\epsilon(p)$ can in principle be solved for arbitrary adiabatic index
$\Gamma(p)$, but here we will spectrally decompose it in terms of a set of
polynomial basis functions,

$\Gamma(p)=\exp\bigg{(}\sum_{k=0}^{n}{\gamma_{k}x^{k}}\bigg{)},$ (24)

where $x=\log(p/p_{0})$ is a dimensionless pressure variable, $p_{0}$ is some
reference pressure, and values of $n$ up to 3 tend to allow for accurate
representations of a variety of EOSs Lindblom (2010).

For given coefficients $\gamma_{k}$ in Eq. (24), the TOV solver available in
LALSuite LIGO Scientific Collaboration (2018) numerically performs the
integrals in Eqs. (23) and (22) to obtain $\epsilon(p)$ and rest mass density
$\rho(p)$, which is inverted to arrive at $p(\rho)$. The reference pressure
$p_{0}$, also called the stitching pressure, is fixed to $5.3716\times
10^{32}$ dyne cm-2; below this pressure, the EOS called SLY Douchin and
Haensel (2001) is stitched on. The coefficients $\gamma_{k}$, $k=0,\ldots,3$
are given uniform priors with ranges $\gamma_{0}\in[0.2,2.0]$,
$\gamma_{1}\in[-1.6,1.7]$, $\gamma_{2}\in[-0.6,0.6]$, and
$\gamma_{3}\in[-0.02,0.02]$. The speed of sound $v_{s}=\sqrt{dp/d\epsilon}$ is
restricted to $v_{s}<1.1\,c$, where the 10% leeway is to allow for imperfect
parameterization. Finally, the adiabatic index is confined to
$\Gamma(p)\in[0.6,4.5]$ Abbott  _et al._ (2018).

###  III.3 Analysis framework

Let $\theta^{a}$ be the components of the parameter vector $\vec{\theta}$,
which in our case will be either the one from Eq. (14) or (15), depending on
whether or not r-modes are taken into account. For GW events with sufficiently
high SNR and assuming stationary, Gaussian noise, the likelihood function for
the signal parameters $\theta^{a}$ will approximately take the form of a
multivariate Gaussian centered on the true values $\hat{\theta}^{a}$ Finn and
Chernoff (1993); Finn (1992). Defining
$\Delta\theta^{a}=\theta^{a}-\hat{\theta}^{a}$, the likelihood becomes

$L(\Delta\theta^{a})=\mathcal{N}e^{-\frac{1}{2}\Gamma_{ab}\Delta\theta^{a}\Delta\theta^{b}},$
(25)

where $\mathcal{N}$ is a normalization factor, and sums over $a$ and $b$ are
implied. The _Fisher matrix_ is given by

$\Gamma_{ab}=\left\langle\frac{\partial
h}{\partial\theta^{a}}\bigg{|}\frac{\partial
h}{\partial\theta^{b}}\right\rangle,$ (26)

where the noise-weighted inner product $\langle\,\cdot\,|\,\cdot\,\rangle$ is
defined as

$\langle
h|g\rangle=4\Re\int_{f_{\text{low}}}^{f_{\text{high}}}\frac{\tilde{h}^{*}(f)\tilde{g}(f)}{S_{n}(f)}df,$
(27)

with $S_{n}(f)$ the one-sided noise power spectral density (PSD).

Strictly speaking, Eq. (25) with (26) pertains to signals $\tilde{h}(f)$ as
seen in a _single_ detector, while the baseline design for Einstein Telescope
assumes three V-shaped detectors arranged in a triangle Punturo  _et al._
(2010); Hild  _et al._ (2011). In principle one would then have to project the
gravitational wave polarizations $\tilde{h}_{+}(f)$, $\tilde{h}_{\times}(f)$
onto each of the three detectors using the appropriate antenna pattern
functions, construct a separate Fisher matrix for each detector, and take the
sum of these to obtain the final Fisher matrix. In this work we will consider
sources from a catalog constructed as in Samajdar  _et al._ (2021), with SNRs
computed for a triangular ET. However, since here we will mainly be interested
in information coming from the phasing, for our purposes it will suffice to
compute a single Fisher matrix with the ET-D PSD Hild  _et al._ (2011) for an
L-shaped detector, and with the SNRs set to the values obtained from the
triangular ET. In Eq. (27) we take the lower frequency cut-off $f_{\rm low}$
to be 5 Hz. For $f_{\rm high}$ we choose the nominal frequency of the
innermosts stable circular orbit (ISCO): $f_{\rm high}=1/(6^{3/2}\pi(1+z)M)$.
We note that for stiff EOSs, the two neutron stars might touch before ISCO is
reached Agathos  _et al._ (2015), but the effect of this for the purposes of
EOS measurements is marginal Wade  _et al._ (2014).

Since we are interested in how well one can measure the EOS with BNS inspiral
signals, we need to choose a “true” EOS, which we take to be FPS Friedman and
Pandharipande (1981), with
$(\gamma_{0},\gamma_{1},\gamma_{2},\gamma_{3})=(1.1561,-0.0468,0.0081,-0.0010)$.
Given a choice of masses $m_{1}$, $m_{2}$, we let the true values of the tidal
deformabilities be $\Lambda_{1}=\Lambda_{\rm FPS}(m_{1})$,
$\Lambda_{2}=\Lambda_{\rm FPS}(m_{2})$, where $\Lambda_{\rm FPS}(m)$ is the
dependence set by the given EOS.

The main question in this work is how well ET will be able to determine the
EOS, i.e., with what accuracy the EOS parameters
$\vec{E}\equiv(\gamma_{0},\gamma_{1},\gamma_{2},\gamma_{3})$ will be measured.
To this end, consider the posterior probability density function
$p(\vec{E},\vec{\theta}^{\prime}|d,\mathcal{H},\mathcal{I})$, where
$\vec{\theta}^{\prime}$ denotes the parameters in (14) or (15) _except_ for
$\tilde{\Lambda}$ and $\delta\tilde{\Lambda}$, which can be calculated from
$\vec{E}$ together with $m_{1}$, $m_{2}$; other than that, $d$ denotes the
data for a given signal, $\mathcal{H}$ a waveform model, and $\mathcal{I}$ any
background information we may possess. Bayes’ theorem tells us that

$p(\vec{E},\vec{\theta^{\prime}}|d,\mathcal{H},\mathcal{I})=\frac{p(d|\vec{E},\vec{\theta}^{\prime},\mathcal{H},\mathcal{I})\,p(\vec{E},\vec{\theta}^{\prime}|\mathcal{H},\mathcal{I})}{p(d|\mathcal{H},\mathcal{I})},$
(28)

where $p(\vec{E},\vec{\theta}^{\prime}|\mathcal{H},\mathcal{I})$ is the prior
probability density for $\vec{E}$ and $\vec{\theta}^{\prime}$,
$p(d|\vec{E},\theta^{\prime},\mathcal{H},\mathcal{I})$ the likelihood, and
$p(d|\mathcal{H},\mathcal{I})$ the evidence, which is set by the requirement
that the posterior probability density be normalized. Clearly
$p(d|\vec{E},\vec{\theta}^{\prime},\mathcal{H},\mathcal{I})$ is not quite the
same as the likelihood $L(\Delta\theta^{a})$ obtained from the Fisher matrix
in Eq. (25), but given true values $\hat{\theta}^{a}$ for the parameters
entering the waveform, it is possible to relate the two by writing

$\displaystyle L(\Delta\theta^{a})$ $\displaystyle=$ $\displaystyle
p(d|\vec{\theta},\mathcal{H},\mathcal{I})$ (29) $\displaystyle=$
$\displaystyle
p(d|\\{\tilde{\Lambda}_{\vec{E}}(m_{1},m_{2}),\delta\tilde{\Lambda}_{\vec{E}}(m_{1},m_{2}),\vec{\theta}^{\prime}\\},\mathcal{H},\mathcal{I})$
$\displaystyle=$ $\displaystyle
p(d|\\{\vec{E},\vec{\theta}^{\prime}\\},\mathcal{H},\mathcal{I}),$

where in the second line, $\tilde{\Lambda}_{\vec{E}}(m_{1},m_{2})$ and
$\delta\tilde{\Lambda}_{\vec{E}}(m_{1},m_{2})$ are the $\tilde{\Lambda}$ and
$\delta\tilde{\Lambda}$ obtained from the component masses for an EOS with the
given parameters $\vec{E}$, and the likelihood in the last line is the one
appearing in Eq. (28).

The way we will proceed is then as follows. Using the Fisher matrix we compute
the likelihood of Eq. (25), and with the identification of Eq. (29) this is
turned into a likelihood in terms of the EOS parameters $\vec{E}$ and waveform
parameters $\vec{\theta}^{\prime}$. We choose flat priors for all of the
individual parameters, and using the emcee sampler Foreman-Mackey  _et al._
(2013) we obtain the posterior density in Eq. (28). By integrating out the
$\vec{\theta}^{\prime}$, this gives us a posterior density for the $\vec{E}$,
$p(\vec{E}|d,\mathcal{H},\mathcal{I})$.111In principle we could have set
$\tilde{\Lambda}=\tilde{\Lambda}_{\vec{E}}(m_{1},m_{2})$ and
$\delta\tilde{\Lambda}=\delta\tilde{\Lambda}_{\vec{E}}(m_{1},m_{2})$ directly
in the Fisher matrix, and obtained 1-sigma uncertainties on the components of
$\vec{E}$ from the covariance matrix in the usual way Finn and Chernoff
(1993); Finn (1992). However, introducing too many additional parameters can
lead to ill-conditioned Fisher matrices, which will cause problems with
numerical inversion. So far we have focused on a single signal, but given a
catalog of detected signals $d_{1},d_{2},\ldots,d_{N}$ that are considered
independent,222If all neutron stars follow the same EOS (as is implicitly
assumed here), then in reality the $d_{1},d_{2},\ldots,d_{N}$ will not be
completely independent. Thus, our way of obtaining combined results may be
somewhat sub-optimal, but this just means that our conclusions about ET’s
ability to reconstruct the EOS will be on the conservative side. information
from all of them can be combined to obtain Del Pozzo  _et al._ (2013)

$\displaystyle p(\vec{E}|d_{1},d_{2},\ldots,d_{N},\mathcal{H},\mathcal{I})$
$\displaystyle=p(\vec{E}|\mathcal{H},\mathcal{I})^{1-N}\prod_{n=1}^{N}p(\vec{E}|d_{n},\mathcal{H},\mathcal{I}),$
(30)

where in practice the posteriors in the product are Gaussian kernel density
estimates of the ones obtained directly with emcee. Through the TOV solver of
LALSuite LIGO Scientific Collaboration (2018), the above combined posterior
distribution for $\vec{E}$ leads to a distribution over equations of state
$p(\rho)$.

Sources are picked from a catalog constructed as in Samajdar  _et al._ (2021).
For binary neutron stars this assumes uniformly distributed source frame
component masses, where for the primary mass $m_{1}\in[1\,M_{\odot},M_{\rm
max}]$ and for the secondary mass $m_{2}\in[1\,M_{\odot},m_{1}]$; here we take
$M_{\rm max}$ to be the maximum mass supported by our reference EOS, which is
$M_{\rm max}=2.03\,M_{\odot}$. As explained in Sec. III.1, we will assume that
for the highest-SNR sources, distances can be measured with sufficient
accuracy that uncertainties on redshift can be neglected in the conversion
between source frame masses and detector frame masses. The sources were
distributed over comoving distance according to a particular prediction for
the merger rate as a function of redshift; for details we refer to the
original paper Samajdar  _et al._ (2021). ET is likely to see tens of
thousands of BNSs per year, but it is reasonable to expect that most of the
information will come from the loudest sources. Hence we consider the 20
loudest sources in the catalog, which have SNRs between 154 and 368. The
angles between the spin vectors and the orbital angular momentum are taken to
be uniform on the sphere. The rotation angular frequencies of the neutron
stars are taken to be uniform in the intervals $\Omega_{i}\in[0,2\pi\times
45]$ Hz, $i=1,2$ Lorimer (2008). Given these ranges for the $\Omega_{i}$ and
the assumed detector lower frequency cut-off of $f_{\rm low}=5$ Hz, there is a
small chance for the r-mode GW frequency $f_{0}$ in Eq. (6) to be below
$f_{\rm low}$ for one or both neutron stars in a binary, but this will not be
the case for any of the simulated sources considered here. We will perform our
analyses with r-modes included, in which case the free parameters are the ones
in Eq. (15), and without r-modes, in which case the free parameters are the
ones in (14), with the values for $\chi_{1z}$ and $\chi_{2z}$ set according to
Eq. (18).

##  IV Results

We are now ready to assess the accuracy with which ET will be able to
reconstruct the EOS of dense nuclear matter. However, in order to check the
validity of our method, we first compare results for a simulated GW170817-like
source seen in Advanced LIGO with those that were obtained in reality for
GW170817. Then we turn to ET and perform an analysis on the 20 loud sources
mentioned above, with or without r-modes included.

###  IV.1 Analysis of a GW170817-like signal

Let us consider a signal with properties similar to that of GW170817, the
binary neutron star signal discovered with LIGO and Virgo in 2017 Abbott  _et
al._ (2017a). Looking at maximum-likelihood parameter values obtained in
Abbott  _et al._ (2019), we take $(m_{1},m_{2})=(1.44,1.27)\,M_{\odot}$,
$\chi_{1z}=0.022$, and $\chi_{2z}=0.0081$. The SNR is set to 32.4, and for the
Fisher matrix we take the PSD to be the one for LIGO Livingston at the time of
the detection, setting $f_{\rm low}=23$ Hz and $f_{\rm high}=1625.36$ Hz. The
effects of r-modes are not included, since they would have had no impact in
this case. Given that we do not know what is the true equation of state, we
consider the same EOS as in the rest of this paper (namely FPS), which for the
given masses yields $(\Lambda_{1},\Lambda_{2})=(387.0,872.2)$.

With the formalism described in the previous section, the above parameters
lead to a “measurement” of $p(\rho)$. Fig. 2 shows the underlying $p(\rho)$
together with a 90% credible region. This is compared with the 90% credible
region for $p(\rho)$ that was actually obtained for this event Abbott  _et
al._ (2018). The qualitative similarity lends confidence to our methodology.

Figure 2: Shown is pressure versus density for a GW170817-like signal whose
EOS is FPS (the red curve), with a 90% credible region (orange). This is
compared with the 90% credible region that was actually obtained for GW170817
(black curves) Abbott  _et al._ (2018). Vertical lines indicate a few
multiples of the nuclear saturation density, $\rho_{\rm sat}$.

###  IV.2 EOS reconstruction with Einstein Telescope

Next we turn to ET. Fig. 3 shows the EOS recovery with only the loudest
source, the 5 loudest sources, and the 20 loudest sources in our simulated
catalog, with r-modes included in the way that was explained in the previous
section. An improvement in measurement accuracy with increasing number of
sources is clearly in evidence. We also note the narrowing of the 90% credible
region near twice the nuclear saturation density, $\rho=2\rho_{\rm sat}$, this
being the approximate average density for most of the neutron stars in our
BNSs, given the EOS we picked. At $2\rho_{\rm sat}$ and for 20 sources, the
width of the 90% credible interval for the pressure is $4.77\times
10^{33}\,\mbox{dyn}/\mbox{cm}^{2}$, to be compared with $8.33\times
10^{34}\,\mbox{dyn}/\mbox{cm}^{2}$ for our simulated GW170817 as seen by
Advanced LIGO. At other densities, ET improves somewhat less on the results
for GW170817; the widths $\Delta P_{90\%}(\rho)$ of the 90% credible intervals
for pressure at a few different values for density are shown in Table 1, for
different numbers of sources in ET, with and without r-modes. Note how at
$2\rho_{\rm sat}$, a single loud source in ET improves the pressure estimation
by a factor of $\sim 3$ over GW170817, but when combining information from 20
sources the improvement is by a factor of $\sim 17$ with r-modes included (and
a factor of $\sim 11$ without r-modes). Still at $2\rho_{\rm sat}$, the gain
from r-modes reaches $\sim 50\%$.

Sources | $\Delta P_{90\%}(\rho_{\rm sat})\,\left[\mbox{dyn}/\mbox{cm}^{2}\right]$ | $\Delta P_{90\%}(2\rho_{\rm sat})\,\left[\mbox{dyn}/\mbox{cm}^{2}\right]$ | $\Delta P_{90\%}(6\rho_{\rm sat})\,\left[\mbox{dyn}/\mbox{cm}^{2}\right]$
---|---|---|---
GW170817-like | $7.08\times 10^{33}$ | $8.33\times 10^{34}$ | $2.07\times 10^{36}$
ET, 1 source |  $3.47\times 10^{33}$ ($3.67\times 10^{33}$) |  $2.62\times 10^{34}$ ($3.33\times 10^{34}$) |  $1.50\times 10^{36}$ ($1.85\times 10^{36}$)
ET, 5 sources |  $2.16\times 10^{33}$ ($2.97\times 10^{33}$) |  $1.16\times 10^{34}$ ($1.79\times 10^{34}$) |  $8.59\times 10^{35}$ ($1.11\times 10^{36}$)
ET, 10 sources |  $8.50\times 10^{32}$ ($1.71\times 10^{33}$) |  $5.94\times 10^{33}$ ($9.37\times 10^{33}$) |  $5.14\times 10^{35}$ ($7.33\times 10^{35}$)
ET, 20 sources |  $6.41\times 10^{32}$ ($1.29\times 10^{33}$) |  $4.81\times 10^{33}$ ($7.54\times 10^{33}$) |  $3.92\times 10^{35}$ ($5.77\times 10^{35}$)
Table 1: The widths of the 90% credible intervals for the pressure at
different densities, for our GW170817-like analysis, and for ET with different
numbers of sources. In the case of ET, the numbers in brackets are without
r-modes. Figure 3: 90% credible regions for pressure versus density from the
loudest source, the 5 loudest sources, and the 20 loudest sources in ET, in
the case where r-modes are included.

Thus, more advanced GW observatories will dramatically improve our knowledge
of the EOS, both through increased sensitivity and by seeing a larger number
of sources. The largest improvement happens near densities that actually occur
in neutron stars. In Fig. 4 we show the accuracy in the measurement of
pressure at $\rho=2\rho_{\rm sat}$, with and without r-modes. The benefit of
seeing r-modes is clearly in evidence. Note how in both cases, most of the
measurement accuracy comes from combining information from the $\sim 10$
loudest sources.

Figure 4: Left panel: The improvement in the measurement of pressure at twice
the nuclear saturation density; shown is the evolution of 90% credible
intervals when going from 1 source to 20 sources, where the blue includes
r-modes while the orange does not. From left to right, sources are being added
in order of decreasing SNR. Right: The individual and combined probability
distributions for pressure.

##  V Summary and conclusions

We have investigated how well one will be able to determine the equation of
state of neutron stars with Einstein Telescope, given that this observatory
will detect tens of BNS inspirals per year for which the SNR will be in the
hundreds. In doing so we have taken into account the effect of resonant
r-modes, which provides additional information about the EOS. The latter were
treated fully relativistically, both in terms of the resonance frequency and
the r-mode overlap; for the former we introduced a new universal relation
linking it to the neutron star tidal deformability, and for the latter we
utilized the recent treatment in Gupta  _et al._ (2021). A reference EOS was
chosen, and general EOSs were represented in terms of the so-called spectral
parameterization. Simulations were performed based on the Fisher approximation
of the likelihood; when identifying the tidal deformabilities $\Lambda_{i}$
with the ones obtained from an EOS with _a priori_ unknown parameters
$\vec{E}$, sampling the likelihood leads to a PDF for $\vec{E}$. From this one
can reconstruct the EOS in terms of pressure as a function of density. We
tested our formalism on a simulated GW170817-like source, and found the
accuracy of our $p(\rho)$ reconstruction to be similar to what was obtained in
reality. We then turned to ET, and focused on the 20 loudest sources in a
simulated “catalog” with realistic assumptions for the merger rate as a
function of redshift.

The estimates arrived at in this paper are necessarily somewhat crude, due to
the approximations made, and especially the limitations of the Fisher matrix
formalism. We only aimed to give a rough sense of how accurately ET would be
able to pin down the EOS; however, the comparison we made for GW170817 between
results from our formalism and actual EOS measurements lends confidence to the
reliability of our predictions.

As expected, ET will dramatically improve our knowledge of the EOS compared
with what was gleaned from GW170817. In terms of $p(\rho)$, the main
improvement comes at densities around twice the nuclear saturation density,
since this is the typical average density of a neutron star. At that density,
ET with 20 loud sources will improve pressure measurements by a factor $\sim
17$ over GW170817. We found that most of the information will come from the
$\sim 10$ loudest sources.

We also saw that the inclusion of r-modes leads to an improvement in our
ability to measure the EOS, at the 50% level in terms of pressure at twice
saturation density. This may seem small compared to what is suggested by Ma
_et al._ (2020) for a single source, but in that work, an optimistic scenario
was assumed where the source had an extremely large SNR (of $\sim 1500$), and
values for the angles $\psi_{1}$, $\psi_{2}$ between the spins and the orbital
angular momentum were chosen so as to nearly maximize the effect of r-modes on
the gravitational waveform. By contrast, our sources had SNRs between 154 and
368, and randomly chosen angles $\psi_{i}$.

In this work we only considered information about the EOS coming from the
inspiral part of the signal. However, ET will also have access to the merger
and post-merger Faber and Rasio (2012), which will lead to further
improvements in the measurement of the EOS. Assessing the effect of the latter
is left for future work.

###### Acknowledgements.

P.K.G., A.P., P.T.H.P., G.K., J.J., and C.V.D.B. are supported by the research
programme of the Netherlands Organisation for Scientific Research (NWO). We
are grateful to Tim Dietrich, Tanja Hinderer, and Jan Steinhoff for useful
discussions at an early stage of the work. This research has made use of data,
software and/or web tools obtained from the Gravitational Wave Open Science
Center (https://www.gw-openscience.org), a service of LIGO Laboratory, the
LIGO Scientific Collaboration and the Virgo Collaboration. LIGO is funded by
the U.S. National Science Foundation. Virgo is funded by the French Centre
National de Recherche Scientifique (CNRS), the Italian Istituto Nazionale
della Fisica Nucleare (INFN) and the Dutch Nikhef, with contributions by
Polish and Hungarian institutes.

## References

  * Aasi  _et al._ (2015) J. Aasi _et al._ (LIGO Scientific), Class. Quant. Grav. 32, 074001 (2015), arXiv:1411.4547 [gr-qc] .
  * Acernese  _et al._ (2015) F. Acernese _et al._ (VIRGO), Class. Quant. Grav. 32, 024001 (2015), arXiv:1408.3978 [gr-qc] .
  * Abbott  _et al._ (2017a) B. P. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. Lett. 119, 161101 (2017a), arXiv:1710.05832 [gr-qc] .
  * Abbott  _et al._ (2020a) B. P. Abbott _et al._ (LIGO Scientific, Virgo), (2020a), arXiv:2001.01761 [astro-ph.HE] .
  * Soares-Santos  _et al._ (2017) M. Soares-Santos _et al._ (DES, Dark Energy Camera GW-EM), Astrophys. J. Lett. 848, L16 (2017), arXiv:1710.05459 [astro-ph.HE] .
  * Cowperthwaite  _et al._ (2017) P. S. Cowperthwaite _et al._ , Astrophys. J. Lett. 848, L17 (2017), arXiv:1710.05840 [astro-ph.HE] .
  * Abbott  _et al._ (2017b) B. P. Abbott _et al._ (LIGO Scientific, Virgo, Fermi-GBM, INTEGRAL), Astrophys. J. Lett. 848, L13 (2017b), arXiv:1710.05834 [astro-ph.HE] .
  * Dietrich  _et al._ (2021) T. Dietrich, T. Hinderer, and A. Samajdar, Gen. Rel. Grav. 53, 27 (2021), arXiv:2004.02527 [gr-qc] .
  * Punturo  _et al._ (2010) M. Punturo _et al._ , Class. Quant. Grav. 27, 084007 (2010).
  * Hild  _et al._ (2011) S. Hild _et al._ , Class. Quant. Grav. 28, 094013 (2011), arXiv:1012.0908 [gr-qc] .
  * Abbott  _et al._ (2017c) B. P. Abbott _et al._ (LIGO Scientific), Class. Quant. Grav. 34, 044001 (2017c), arXiv:1607.08697 [astro-ph.IM] .
  * Reitze  _et al._ (2019) D. Reitze _et al._ , Bull. Am. Astron. Soc. 51, 035 (2019), arXiv:1907.04833 [astro-ph.IM] .
  * Samajdar  _et al._ (2021) A. Samajdar, J. Janquart, C. Van Den Broeck, and T. Dietrich, Phys. Rev. D 104, 044003 (2021), arXiv:2102.07544 [gr-qc] .
  * Maggiore  _et al._ (2020) M. Maggiore _et al._ , JCAP 03, 050 (2020), arXiv:1912.02622 [astro-ph.CO] .
  * Evans  _et al._ (2021) M. Evans _et al._ , (2021), arXiv:2109.09882 [astro-ph.IM] .
  * Flanagan and Hinderer (2008) E. E. Flanagan and T. Hinderer, Phys. Rev. D 77, 021502 (2008), arXiv:0709.1915 [astro-ph] .
  * Hinderer (2008) T. Hinderer, Astrophys. J. 677, 1216 (2008), arXiv:0711.2420 [astro-ph] .
  * Del Pozzo  _et al._ (2013) W. Del Pozzo, T. G. F. Li, M. Agathos, C. Van Den Broeck, and S. Vitale, Phys. Rev. Lett. 111, 071101 (2013), arXiv:1307.8338 [gr-qc] .
  * Agathos  _et al._ (2015) M. Agathos, J. Meidam, W. Del Pozzo, T. G. F. Li, M. Tompitak, J. Veitch, S. Vitale, and C. Van Den Broeck, Phys. Rev. D 92, 023012 (2015), arXiv:1503.05405 [gr-qc] .
  * Lackey and Wade (2015) B. D. Lackey and L. Wade, Phys. Rev. D 91, 043002 (2015), arXiv:1410.8866 [gr-qc] .
  * Abbott  _et al._ (2018) B. P. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. Lett. 121, 161101 (2018), arXiv:1805.11581 [gr-qc] .
  * Abbott  _et al._ (2020b) B. P. Abbott _et al._ (LIGO Scientific, Virgo), Class. Quant. Grav. 37, 045006 (2020b), arXiv:1908.01012 [gr-qc] .
  * Flanagan and Racine (2007) E. E. Flanagan and E. Racine, Phys. Rev. D 75, 044001 (2007), arXiv:gr-qc/0601029 .
  * Poisson (2020) E. Poisson, Phys. Rev. D 101, 104028 (2020), arXiv:2003.10427 [gr-qc] .
  * Ma  _et al._ (2020) S. Ma, H. Yu, and Y. Chen, (2020), arXiv:2010.03066 [gr-qc] .
  * Idrisy  _et al._ (2015) A. Idrisy, B. J. Owen, and D. I. Jones, Phys. Rev. D 91, 024001 (2015), arXiv:1410.7360 [gr-qc] .
  * Gupta  _et al._ (2021) P. K. Gupta, J. Steinhoff, and T. Hinderer, Phys. Rev. Res. 3, 013147 (2021), arXiv:2011.03508 [gr-qc] .
  * Jiménez Forteza  _et al._ (2018) X. Jiménez Forteza, T. Abdelsalhin, P. Pani, and L. Gualtieri, Phys. Rev. D 98, 124014 (2018), arXiv:1807.08016 [gr-qc] .
  * LIGO Scientific Collaboration (2018) LIGO Scientific Collaboration, “LIGO Algorithm Library - LALSuite,” free software (GPL) (2018).
  * Sathyaprakash and Dhurandhar (1991) B. S. Sathyaprakash and S. V. Dhurandhar, Phys. Rev. D44, 3819 (1991).
  * Blanchet (2002) L. Blanchet, Living Rev. Rel. 5, 3 (2002), arXiv:gr-qc/0202016 .
  * Vines  _et al._ (2011) J. Vines, E. E. Flanagan, and T. Hinderer, Phys. Rev. D83, 084051 (2011), arXiv:1101.1673 [gr-qc] .
  * Kidder  _et al._ (1993) L. E. Kidder, C. M. Will, and A. G. Wiseman, Phys. Rev. D 47, R4183 (1993), arXiv:gr-qc/9211025 .
  * Kidder (1995) L. E. Kidder, Phys. Rev. D 52, 821 (1995), arXiv:gr-qc/9506022 .
  * Apostolatos  _et al._ (1994) T. A. Apostolatos, C. Cutler, G. J. Sussman, and K. S. Thorne, Phys. Rev. D 49, 6274 (1994).
  * Zhao and Wen (2018) W. Zhao and L. Wen, Phys. Rev. D 97, 064031 (2018), arXiv:1710.05325 [astro-ph.CO] .
  * Ade  _et al._ (2014) P. Ade _et al._ (Planck), Astron. Astrophys. 571, A16 (2014), arXiv:1303.5076 [astro-ph.CO] .
  * Wade  _et al._ (2014) L. Wade, J. D. E. Creighton, E. Ochsner, B. D. Lackey, B. F. Farr, T. B. Littenberg, and V. Raymond, Phys. Rev. D 89, 103012 (2014), arXiv:1402.5156 [gr-qc] .
  * Yagi and Yunes (2013) K. Yagi and N. Yunes, Science 341, 365 (2013), arXiv:1302.4499 [gr-qc] .
  * Koliogiannis and Moustakidis (2019) P. S. Koliogiannis and C. C. Moustakidis, Astrophys. Space Sci. 364, 52 (2019), arXiv:1806.09999 [nucl-th] .
  * Lindblom (2010) L. Lindblom, Phys. Rev. D 82, 103011 (2010), arXiv:1009.0738 [astro-ph.HE] .
  * Douchin and Haensel (2001) F. Douchin and P. Haensel, Astron. Astrophys. 380, 151 (2001), arXiv:astro-ph/0111092 .
  * Finn and Chernoff (1993) L. S. Finn and D. F. Chernoff, Phys. Rev. D 47, 2198 (1993), arXiv:gr-qc/9301003 .
  * Finn (1992) L. S. Finn, Phys. Rev. D 46, 5236 (1992).
  * Friedman and Pandharipande (1981) B. Friedman and V. R. Pandharipande, Nucl. Phys. A 361, 502 (1981).
  * Foreman-Mackey  _et al._ (2013) D. Foreman-Mackey, D. W. Hogg, D. Lang, and J. Goodman, Publications of the Astronomical Society of the Pacific 125, 306–312 (2013).
  * Lorimer (2008) D. R. Lorimer, Living Rev. Rel. 11, 8 (2008), arXiv:0811.0762 [astro-ph] .
  * Abbott  _et al._ (2019) B. P. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. X 9, 011001 (2019), arXiv:1805.11579 [gr-qc] .
  * Faber and Rasio (2012) J. A. Faber and F. A. Rasio, Living Rev. Rel. 15, 8 (2012), arXiv:1204.3858 [gr-qc] .


# Sterile Neutrino Sensitivity with Wrong-Sign Muon Appearance at $\nu$STORM

C.D. Tunnell

IDS-NF-035

Figure 1: A diagram of the proposed $\nu$STORM facility.

Neutrinos from STORed Muons111The facility was previously called the Very Low
Energy Neutrino Factory (VLENF). ($\nu$STORM) is a proposed experiment that
uses 3.8 GeV/c muon decay to produce a well-understood beam of electron and
muon neutrinos that can be used for short baseline physics (Fig. 1). A
magnetized far detector allows for the wrong-sign muon appearance physics of
$\nu_{e}\to\nu_{\mu}$ and provides more sensitivity to sterile neutrinos than
other proposals (See comparisons in [1]) because of the relative ease to which
muon tracks can be identified. Other physics such as $\nu_{e}$ and $\nu_{\mu}$
cross section measurements are possible. For further details, see the Letter
of Intent [16].

An explanation of the $\nu$STORM appearance analysis will follow. This work is
a continuation of the work presented in [20]. For disappearance measurement
work, see [21].

##  1 Short Baseline Oscillations

LEP experiments revealed that there are three light neutrinos that couple to
the $Z$-boson (ie. _active neutrinos_), however, there are theoretical and
experimental motivations [1] for neutrinos without Standard Model interactions
called _sterile_ neutrinos. The (3+1) scenario is the case of three active
neutrinos with an additional heavy sterile neutrino –
$m_{4}>>m_{\text{others}}$ – and only this situation is considered although
the results are generalizable.

The probability $\nu_{e}\to\nu_{\mu}$ depends on the mixing matrix $U$
(Reviewed in [5]). Let $R_{ij}$ be a rotation between the $i$-th and $j$-th
mass eigenstates without a CP violating phase: CP violation cannot be observed
in oscillations with large $\Delta m^{2}$ dominance (See p.g. 273 of [7]). For
$N$ neutrinos, $R_{ij}$ has dimension $N\times N$ and takes the form:

$\displaystyle R_{ij}=\begin{pmatrix}1&\ldots&0&\ldots&0&\ldots&0\\\
\vdots&&\vdots&&\vdots&&\vdots\\\
0&\ldots&\cos\theta_{ij}&\ldots&\sin\theta_{ij}&\ldots&0\\\
\vdots&&\vdots&&\vdots&&\vdots\\\
0&\ldots&-\sin\theta_{ij}&\ldots&\cos\theta_{ij}&\ldots&0\\\
\vdots&&\vdots&&\vdots&&\vdots\\\ 0&\ldots&0&\ldots&0&\ldots&1\\\
\end{pmatrix}.$ (1)

By convention, the three neutrino mixing matrix is
$U_{\text{PMNS}}=R_{23}R_{13}R_{12}$. In the (3+1) model of neutrino
oscillations, extra rotations can be introduced such that the mixing matrix is
$U_{\text{(3+1)}}=R_{34}R_{24}R_{14}U_{\text{PMNS}}$. Given that $\Delta
m^{2}_{41}>>\Delta m^{2}_{31}$, $U_{\text{PMNS}}$ can be approximated by the
identity matrix (_ie._ the “short baseline approximation”) implying
$U_{e4}=\sin(\theta_{14})$ and $U_{\mu 4}=\sin(\theta_{24})\cos(\theta_{14})$.

The oscillation probabilities for appearance and disappearance, respectively,
are:

$\displaystyle\text{P}_{\nu_{e}\to\nu_{\mu}}$ $\displaystyle=$ $\displaystyle
4|U_{e4}|^{2}|U_{\mu 4}|^{2}\sin^{2}\left(\frac{\Delta
m^{2}_{41}L}{4E}\right)$ (2) $\displaystyle=$
$\displaystyle\sin^{2}(2\theta_{e\mu})\sin^{2}\left(\frac{\Delta
m^{2}_{41}L}{4E}\right),$ (3)
$\displaystyle\text{P}_{\nu_{\alpha}\to\nu_{\alpha}}$ $\displaystyle=$
$\displaystyle 1-\left[4|U_{\alpha 4}|^{2}(1-|U_{\alpha
4}|^{2})\right]\sin^{2}\left(\frac{\Delta m^{2}_{41}L}{4E}\right).$ (4)

in this short baseline limit where the definition
$\sin^{2}(2\theta_{e\mu})=4|U_{e4}|^{2}|U_{\mu 4}|^{2}$ has been introduced.

Electron and muon neutrino disappearance measurements will constrain
$|U_{e4}|^{2}$ ([21]) and $|U_{\mu 4}|^{2}$ while the appearance channel
analysis could measure the product $|U_{e4}|^{2}|U_{\mu 4}|^{2}$. Information
about the matrix element $U_{e4}$ also arises from jointly analyzing
$\bar{\nu}_{\mu}$ disappearance and $\nu_{\mu}$ appearance. The remaining
matrix element $U_{\tau 4}$ can be extracted by analyzing NC rates
$|U_{s4}|^{2}=\sum_{e,\mu,\tau}|U_{\alpha 4}|^{2}$, using the other channels
to constrain $|U_{e4}|^{2}$ and $|U_{\mu 4}|^{2}$, and assuming unitarity.

##  2 The Neutrino Flux: $\Phi$

Table 1: Matrix elements for muon decay | $f_{0}(x)$ | $f_{1}(x)$
---|---|---
$\nu_{\mu}$ | $2x^{2}(3-2x)$ | $2x^{2}(1-2x)$
$\nu_{e}$ | $12x^{2}(1-x)$ | $12x^{2}(1-x)$

Muon-decay beams contrast pion-decay beams because the beam characteristics
and production mechanisms are well-known. The neutrino flux arises from the
Electroweak decay of $\mu\to\nu_{\mu}\bar{\nu}_{e}e$ and it is sufficient to
compute matrix elements at tree level. The neutrino spectrum for a
$\mu^{\pm}\to e^{\pm}+\nu_{e}(\bar{\nu}_{e})+\bar{\nu_{\mu}}(\nu_{\mu})$ decay
in the rest frame of the muon follows:

$\displaystyle\frac{\operatorname{d}n}{\operatorname{d}x\operatorname{d}\Omega}=\frac{1}{4\pi}\left[f_{0}(x)\mp\mathcal{P}f_{1}(x)\cos\theta\right]$
(5)

where $x=2E^{\text{c.o.m.}}_{\nu}/m_{\mu}\in[0,1]$ is the scaled neutrino
energy in the rest frame, $\Omega$ is the solid angle in the rest frame,
$f_{0}(x)$ and $f_{1}(x)$ are muon decay parameters, and $\mathcal{P}$ is the
polarization. Electron and neutrino masses are negligible for this process and
ignored, hence the inclusive range for values of $x$. These muon decay
paramters can be computed to leading order with Electroweak theory (See, for
example, chapter 6 of Ref. [18]) and are neutrino flavor dependent (See Table
1).

Figure 2: The flux of $\nu_{e}$ and $\bar{\nu}_{\mu}$ for a 3.8GeV/c muon-
decay without oscillations at 2000 meters. No smearing due to accelerator
effects has been performed.

The polarization $\mathcal{P}$ is set to zero, similar to other studies, and
has been shown to average to zero due to Thomas Precession. Boosting the
neutrino distributions into the lab frame leads to:

$\displaystyle\frac{\operatorname{d}^{2}N_{\mu}}{\operatorname{d}y\operatorname{d}A}$
$\displaystyle=$ $\displaystyle\frac{4n_{\mu}}{\pi
L^{2}m^{6}_{\mu}}E^{4}_{\mu}y^{2}(1-\beta\cos\phi)\left[3m^{2}_{\mu}-4E^{2}_{\mu}y(1-\beta\cos\phi)\right]$
(6)
$\displaystyle\frac{\operatorname{d}^{2}N_{e}}{\operatorname{d}y\operatorname{d}A}$
$\displaystyle=$ $\displaystyle\frac{24n_{\mu}}{\pi
L^{2}m^{6}_{\mu}}E^{4}_{\mu}y^{2}(1-\beta\cos\phi)\left[m^{2}_{\mu}-2E^{2}_{\mu}y(1-\beta\cos\phi)\right]$
(7)

where $y=E_{\nu}/E_{\mu}$ is the scaled neutrino energy in the lab frame,
$\beta=\sqrt{1-m^{2}_{\mu}/E^{2}_{\mu}}$, $A$ is an area, and $n_{\mu}$ is the
number of muons. These neutrino distributions (Fig. 2) are for a point source
so they are not directly applicable to the decay straight of $\nu$STORM.

The number of muons assumed is $1.8\times 10^{18}$ and is based on $10^{21}$
protons on target (POT) at 60 GeV/c. It corresponds to roughly 5 years of
running with a 100 kW target station. The number of useful muon decays is
motivated in [16].

Figure 3: The unoscillated flux of $\nu_{e}$ and $\bar{\nu}_{\mu}$ for a
$(3.8\pm 0.38)\text{ GeV/c}$ muon-decay at 2000 meters. Accelerator effects
are included; see the text for details. Figure 4: The flux at the far detector
for a $(3.8\pm 0.38)\text{ GeV/c}$ muon for initial $\nu_{e}$ states including
integration over the beam envelope and detector volume. Final states include
$\nu_{e}$ without oscillations and both $\nu_{e}$ and $\bar{\nu}_{\mu}$ with
best fit short baseline oscillations. The normalization is $10^{21}$ POT.

When computing the flux for $\nu$STORM, the far detector approximation of a
point-source accelerator and detector no longer is applicable since the size
of the detector and accelerator straight (150 meters) are comparable to the
baseline of 2000 meters. The neutrino fluxes are computed by integrating over
the decay straight, transverse beam phase space, and detector volume. The beam
occupies a 6D phase space ($x$, $y$, $z$, $p_{x}$, $p_{y}$, $p_{z}$) and the
detector has a $5\text{ m}\times 5\text{ m}$ cross section with the depth set
by the desired fiducial mass of 1.3 kt. Both transverse 2D phase spaces are
represented by the Twiss parameters $\alpha=0$ and $\beta=40\text{ m}$ where
the $1\sigma$ Gaussian geometric emittance is assumed to be $2.1\text{ mm}$.
The spread in, for example, $x$ is $\sigma_{x}=\sqrt{\beta\epsilon}$ and the
angular divergence in $x$ is $\sigma_{x^{\prime}}=\sqrt{\epsilon/\beta}$. The
longitudinal phase space ($z$ and $p_{z}$) is described by assuming a uniform
distribution in $z\in[0,150\text{ m}]$ – accurate to 0.5% – and
$p_{z}\in[3.8\pm 0.38\text{ GeV/c}]$.

The flux is computed by Monte Carlo (MC) integration: random points are chosen
within the beam phase space and within the detector volume to determine the
expected flux. This integration introduces a new computational requirement:
the baseline is a variable that affects both the oscillation probability
($L/E$) and the flux ($L^{-2}$ geometric factor). The GLoBES software (version
3.1.10) [12, 11] that is used for neutrino factory phenomenology treats these
as separable problems and was modified to compute this flux (and later the
event rates and sensitivities). Specifically, GLoBES is modified such that
both the flux and oscillation probability are computed in the _oscillation
probability engine_. The code for the analysis is available [19] under the GPL
license [9].

The resulting flux after the integration (Fig. 3) is corrected for accelerator
effects. The corrections are small for far detector physics (Compare to Fig.
2) but are important for near detector physics where the baseline is smaller
than the decay straight.

##  3 The Oscillation Probability: (Prob.)

Figure 5: The oscillation probability for the “golden channel”
$\nu_{e}\to\nu_{\mu}$ from Eq. 2 using the (3+1) oscillation parameters in
TABLE 2. A baseline of 2000 meters is assumed.

This section will discuss how sterile oscillation phenomenology relates to
conducting the proposed experiment. For instance, for a point-source baseline
of 2000 meters, it is possible to determine the oscillation probability (Fig.
5) using Eq. 2 for any combination of $L$ and $E$.

Table 2: Best-fit oscillation parameters for the (3+1) sterile neutrino scenario using combined MB $\bar{\nu}$ and LSND $\bar{\nu}$ data [8]. Parameter | Value
---|---
$\Delta m^{2}_{41}$ [$\text{eV}^{2}$] | 0.89
$|U_{e4}|^{2}$ | 0.025
$|U_{\mu 4}|^{2}$ | 0.023
Table 3: Values for $3\times 3$ oscillations used. $\sin^{2}\theta_{12}=0.319$
---
$\sin^{2}\theta_{23}=0.462$
$\sin^{2}\theta_{13}=0.010$
$\Delta m^{2}_{21}=7.59\times 10^{-5}\text{ eV}^{2}$
$\Delta m^{2}_{31}=2.46\times 10^{-3}\text{ eV}^{2}$

The best fit parameters for the “short baseline anomaly” and $3\times 3$
mixing (_i.e._ $\sin^{2}(2\theta_{13})$, $\Delta m^{2}_{12}$, etc.) are used
throughout the analysis. The best fit parameters for the LSND anomaly come
from [8] (See TABLE 2) and agree with those published by the LSND
collaboration [2]. For completeness, oscillations between known mass
eigenstates are included despite not influencing the sensitivity: the
correction is order $10^{-5}$. The best fit data from [10] is used to specify
standard $3\times 3$ oscillations. Without loss of generality, normal
hierarchy is assumed and the values of known $3\times 3$ mixing can be seen in
Table 3. Errors associated with these quantities are ignored.

Computationally, the SNU (version 1.1) add-on [13, 14] has been used to extend
computations in GLoBES to $4\times 4$ mixing matrices.

##  4 Cross section: $\sigma$

(a) CC

(b) NC

Figure 6: Neutrino cross sections per nucleon.

Cross sections are required for each neutrino flavor ($\nu_{\mu}$,
$\bar{\nu}_{\mu}$, $\nu_{e}$, $\bar{\nu}_{e}$) and each interaction type (CC
or NC). The nucleon cross sections (Fig. 6) are calculated in [15] and [17]
for the low energy and high energies, respectively. NC cross sections are
flavor independent. The CC cross sections are approximately flavor
independent: _Fermi’s Second Golden Rule_ results in the same matrix elements
and, at these energies, the phase spaces for the final-state electrons and
muons are equal.

The total cross section requires knowing the number of nucleons in addition to
the nucleon cross section. The fiducial mass of 1.3 kt determines the number
of nucleons via Avogadro’s number.

##  5 Interaction rates: $N_{\text{int.}}$

(a) Appearance with stored $\mu^{+}$

(b) Appearance with stored $\mu^{-}$

(c) Disappearance with stored $\mu^{+}$

(d) Disappearance with stored $\mu^{-}$

Figure 7: True channel rate energy distributions assuming the LSND anomaly
best fit values. The transitions $\nu_{e}\to\nu_{\mu}$,
$\bar{\nu}_{e}\to\bar{\nu}_{\mu}$, $\bar{\nu}_{\mu}\to\bar{\nu}_{\mu}$, and
$\nu_{\mu}\to\nu_{\mu}$ are shown.

The number of neutrino interactions is computed which does not require
assumptions about the detector. The interaction rates can be computed by
$N_{\text{int.}}=\Phi\times\text{(Prob.)}\times\sigma$, for flux $\Phi$,
oscillation probability $(\text{Prob.})$, and cross section $\sigma$, where
all of these quantities have been computed in the previous sections.

Using the LSND anomaly best fit (TABLE 2) as a example for a sterile neutrino
signal, the event rates for $\mu^{+}$ and $\mu^{-}$ decays are shown in TABLE
4. Various deductions can be made about these event rates and their
statistical significance. With either stored $\mu^{+}$s or stored $\mu^{-}$s,
the statistical significance of all channels is greater than $10\sigma$.
Combining the NC channels together results in a statistical significance of
$20\sigma$ and $17\sigma$ for stored $\mu^{+}$ and $\mu^{-}$, respectively.
There are no known physics backgrounds to neither $\nu_{e}\to\nu_{m}u$ CC nor
$\bar{nu}_{e}\to\bar{\nu}_{m}u$ CC interactions except to negligible solar-
term oscillations, so the backgrounds will arise from how well the detector
can differentiate these interactions.

The number of events can also be determined as a function of energy since the
evolution of $\rho$, $\sigma$, and $(\text{Prob.})$ as a function of energy is
known. These distributions are shown in Fig. 7.

There are numerous channels with reach into the sterile neutrino parameter
space. Most other experiments have one channel to explore (See [1] for list of
experiments), whereas in the best case $\nu$STORM allows for 10 signals and in
the worst case 6 (i.e. combine $\nu_{e}\to\nu_{e}$ CC and all NC channels).

Table 4: Truth event rates for $10^{21}$ POT for the no oscillations and short
baseline oscillations described by TABLE 2. The statistical significances are
computed. The combined statistical significance of NC events are 20 and 17 for
stored $\mu^{+}$ and $\mu^{-}$, respectively. There are no physics backgrounds
to $\nu_{e}\to\nu_{m}u$ CC interactions.

Channel | $N_{\textrm{osc.}}$ | $N_{\textrm{null}}$ | Diff. | $(N_{\textrm{osc.}}-N_{\textrm{null}})/\sqrt{N_{\textrm{null}}}$
---|---|---|---|---
$\bar{\nu}_{e}\to\bar{\nu}_{\mu}$ CC | 117 | 0 | $\infty$ | $\infty$
$\bar{\nu}_{e}\to\bar{\nu}_{e}$ NC | 30511 | 32481 | -6.1% | -10.9
$\nu_{\mu}\to\nu_{\mu}$ NC | 66037 | 69420 | -4.9% | -12.8
$\bar{\nu}_{e}\to\bar{\nu}_{e}$ CC | 77600 | 82589 | -6.0% | -17.4
$\nu_{\mu}\to\nu_{\mu}$ CC | 197284 | 207274 | -4.8% | -21.9
(a) Stored $\mu^{-}$.

Channel | $N_{\textrm{osc.}}$ | $N_{\textrm{null}}$ | Diff. | $(N_{\textrm{osc.}}-N_{\textrm{null}})/\sqrt{N_{\textrm{null}}}$
---|---|---|---|---
$\nu_{e}\to\nu_{\mu}$ CC | 332 | 0 | $\infty$ | $\infty$
$\bar{\nu}_{\mu}\to\bar{\nu}_{\mu}$ NC | 47679 | 50073 | -4.8% | -10.7
$\nu_{e}\to\nu_{e}$ NC | 73941 | 78805 | -6.2% | -17.3
$\bar{\nu}_{\mu}\to\bar{\nu}_{\mu}$ CC | 122322 | 128433 | -4.8% | -17.1
$\nu_{e}\to\nu_{e}$ CC | 216657 | 230766 | -6.1% | -29.4
(b) Stored $\mu^{+}$.

##  6 Event rates after cuts

It must be determined how many of the raw events pass analysis cuts. Similar
analyses have been performed for Neutrino Factories exploring CP violation at
energies ranging from 25 GeV [4] to 5 GeV [6], but never at 3.8 GeV/c.
Preexisting experience and knowledge exists as to fractional background levels
and analysis difficulties; work had to be performed in order to tune the
analysis for this energy range.

The detector performance can be represented by _migration matrices_ (also
known as response matrices or energy smearing matrices) that describe both the
energy resolution and detection efficiency. If events are binned in terms of
true neutrino energy then the migration matrix is needed to transform the
distribution into the space of measured neutrino energies. For example, take
the histogram:

$\vec{h}^{\text{true}}=(N^{\text{true}}_{\text{0.0 - 0.1
\text{GeV}}},N^{\text{true}}_{\text{0.1 - 0.2
\text{GeV}}},\ldots,N^{\text{true}}_{\text{3.9 - 4.0 \text{GeV}}})^{T},$ (8)

where $N^{\text{true}}_{\text{0.0 - 0.1 \text{GeV}}}$ is the number of events
in the bin with ranges 0.0 and 0.1 GeV. The migration matrix $\mathbf{M}$ used
for this analysis is a square matrix and defined such that
$\vec{h}^{\text{measured}}=\mathbf{M}\vec{h}^{\text{true}}$ where
$\vec{h}^{\text{measured}}$ is the expected histogram of reconstructed
quantities in the detector.

With a perfect detector $\mathbf{M}=\text{diag.}(1,1,\ldots,1)$. $\mathbf{M}$
is unitary if and only if it describes only energy smearing. Efficiencies are
included into $\mathbf{M}$ by removing the unitarity constraint.

(a) $\nu_{\mu}$ CC

(b) $\bar{\nu}$ NC

(c) $\bar{\nu}_{\mu}$ CC

(d) $\nu_{e}$ CC

Figure 8: Migration matrices for $\nu_{\mu}$ CC, $\bar{\nu}$ NC,
$\bar{\nu}_{\mu}$ CC, and $\nu_{e}$ CC.

Migration matrices have been computed for $\nu_{\mu}$ CC, $\bar{\nu}_{\mu}$
CC, $\bar{\nu}_{\mu}$ NC, and $\nu_{e}$ CC (See [16]) and can be seen in Fig.
8. The background level of $\nu_{e}$ NC events into the signal window are
negligible compared with $\bar{\nu}_{\mu}$ NC due to the lower energies. These
numbers are derived using a GENIE and Geant4 simulation, described in the
cited text, which are MC method softwares. Statistical fluctuations exist in
the migration matrices due to computational limitations.

Figure 9:  The rule rate as a function of observed energy for the appearance
channel $\nu_{e}\to\nu_{\mu}$. Migration matrices are used for $\nu_{\mu}$ CC,
$\bar{\nu}_{\mu}$ CC, $\bar{\nu}_{\mu}$ NC, and $\nu_{e}$ CC.

##  7 Statistics

It is necessary to determine if the number of events observed after cuts
(_i.e._ rule rates) is statistically significant. The experiment must reject
the null hypothesis when accounting for statistical fluctuations.

The hypothesis $H_{0}$ of no oscillations is the null hypothesis and designate
$H_{1}$ to be the alternate hypothesis. These hypotheses have oscillation
parameters associated with them: let $\theta_{0}=\\{\Delta
m^{2}_{41},\theta_{34},\theta_{24},\theta_{14}\\}$ be the oscillation
parameters associated with $H_{0}$, and similarly $\theta_{1}$ for $H_{1}$.

The _test statistic_ $X$ is a function of the experimental observations and
let $W$ be the space of all possible values of $X$. One can divide $W$ into
two regions: the region $w$ for those possible values of $X$ which would
suggest that the null hypothesis $H_{0}$ is not true and the remaining region
$W-w$.

It is desirable to have a small probability of $X$ – by statistical
fluctuations alone – taking a value in $w$ when $H_{0}$ is true. A level of
significance $\alpha$ can be defined:

$P(X\in w|H_{0})=\alpha$ (9)

where $\alpha$ corresponds to, colloquially, “5$\sigma$” when $\alpha\simeq
2.8\times 10^{-7}$ and “10$\sigma$” when $\alpha\simeq 7.6\times 10^{-24}$.
The number of “$\sigma$” correspond to the $p$-value of having a greater than
$n\sigma$ upward fluctuation of a Gaussian centered at zero. No Gaussian
assumptions are made in this analysis.

The test statistic that will be used for hypothesis testing is the likelihood
ratio test. Let there be $N$ observations $\mathbf{X}=\\{X_{1},...,X_{N}\\}$
and a probability distribution function $f(X_{i}|\theta)$. The likelihood
function is:

$\displaystyle L(\mathbf{X}|\mathbf{\theta})$ $\displaystyle=$
$\displaystyle\prod^{N}_{i=1}f(X_{i}|\mathbf{\theta})$ (10) $\displaystyle=$
$\displaystyle\prod_{i}e^{-\lambda_{i}}\lambda^{X_{i}}_{i}/{X_{i}}!$ (11)

where $\lambda_{i}$ is the expected number of background in the bin with
$X_{i}$ events and is a function of $\theta$. The distribution is Poisson
because the background levels are small. The short baseline parameters
$\theta_{1}$ for $H_{1}$ are free to take any value but the parameters
$\theta_{0}$ are fixed to zero by the null hypothesis requiring no
oscillations. The likelihood ratio test defines a test statistic $\lambda$
such that:

$\lambda=\frac{L(\mathbf{X}|\mathbf{\theta_{0}})}{\max_{\theta_{1}}L(\mathbf{X}|\mathbf{\theta_{1}})}$
(12)

where the denominator is maximized with respect to $\theta_{1}$ while the
numerator remains fixed. Using Eq. 11 leads to:

$\lambda=\prod_{i}e^{-\lambda_{i}+X_{i}}\left(\lambda_{i}/X_{i}\right)^{X_{i}}.$
(13)

The $\chi^{2}$ can be defined as $\chi^{2}=-2\ln\lambda$ (See [3]) which is
preferable to using $\lambda$ because of specifics about how multiplication is
performed by a computer. Using this definition, one finds:

$\displaystyle\chi^{2}=-2\ln\lambda=2\sum_{i}\lambda_{i}-X_{i}+X_{i}\ln\left(\frac{\lambda_{i}}{X_{i}}\right)$
(14)

which has two degrees of freedom since the numerator of Eq. 12 has no degrees
of freedom and the denominator has two degrees of freedom.

A test statistic has been defined that allows for determining if an experiment
is sensitive to various oscillation parameters. The $\chi^{2}$ can be computed
in terms of energy bins, with the appropriate definition of $X_{i}$, allowing
for spectral information to be used when computing sensitivities.

##  8 The Appearance Analysis

The parameters to be explored in the appearance analysis are $\Delta
m^{2}_{41}$ and $\sin^{2}(\theta_{e\mu})$. Contours in the neutrino parameter
space $\Delta m^{2}_{41}$ versus $\sin^{2}(\theta_{e\mu})$ can be used to
compare the sensitivities of various proposed short baseline experiments. A
statistics-only $\chi^{2}$ using spectral information is used (Fig.10).

Care must be taken when defining $\chi^{2}(\Delta
m^{2}_{41},\sin^{2}(\theta_{e\mu}))$ to ensure that it is well-defined. In the
(3+1) scenario, the signal $\nu_{e}\to\nu_{\mu}$ depends on the amplitude
$\sin^{2}(\theta_{e\mu})=4|U_{e4}|^{2}|U_{\mu 4}|^{2}$ and frequency $\Delta
m^{2}_{41}$ (See Eq. 2). If there is an appearance signal, then
$|U_{e4}|^{2}|U_{\mu 4}|^{2}\neq 0$ which implies that both $U_{e4}$ and
$U_{\mu 4}$ are nonzero. There is disappearance of the CC and NC backgrounds
(See Eq. 4) which affects the background estimation in the $\chi^{2}$. This
issue is addressed by not oscillating the backgrounds thus overestimating the
backgrounds.

Figure 10: Sterile sensitivity under the appearance channel
$\nu_{e}\to\nu_{\mu}$. This channel is the CPT of the LSND anomaly
$\bar{\nu}_{\mu}\to\bar{\nu}_{e}$. There is $10\sigma$ sensitivity to the LSND
and MiniBooNe 99% confidence interval [8]. Figure 11:  A baseline optimization
using a total rates statistics-only $\chi^{2}$, a signal efficiency of 0.5,
and background rejection of charge misidentification and NCs at $10^{-3}$ and
$10^{-4}$. Figure 12:  Tuning the NC rejection cut. The NC rejection level is
shown versus the signal efficiency. A charge misidentification background of
$10^{-4}$ is shown to illustrate when NC backgrounds become statistically
significant. A total rates statistics-only $\chi^{2}$ is used. Figure 13:
Tuning the charge misidentification cut. The charge misidentification level is
shown versus the signal efficiency. A NC background of $10^{-4}$ is shown to
illustrate when charge misidentification backgrounds become statistically
significant. A total rates statistics-only $\chi^{2}$ is used. Figure 14:  An
optimization between the detector performance and accelerator performance
using the charge misidentification rates and number of muon decays as the
performance metric. IDR refers to the Interim Design Report [4] detector
performance. FODO refers to the FODO lattice design that gives $1.8\times
10^{18}$ useful muon decays whilst FFAG refers to the FFAG design that gives
$4.68\times 10^{18}$ useful muon decays. Both accelerators assume a front-end
of the main injector at 60 GeV/c.

As the cuts-based detector performance section improves and various cost
optimizations are done, there are numerous parameters that can be tuned to
compensate and conserve the physics that can be done with such a facility. For
example, the optimization of baseline and energy (Fig. 11) allows one to
change the baseline depending on site constraints or modify the energy of the
ring if the accelerator gets too expensive. As the cuts-based detector
performance improves, the various background rejections (Fig. 12 and 13) may
allow for a smaller detector or cheaper target station. The tools have been
developed that allow the important accelerator and detector performance
metrics into cost optimizations.

Figure 11 shows that, for a fixed baseline, increasing the muon energy is
always advantageous. This effect arises because the maximum of the $\nu_{e}$
flux is not at the oscillation maximum but rather at a higher energy. At high
energies the oscillation probability is:

$\displaystyle\text{Pr}[\nu_{e}\to\nu_{\mu}]$ $\displaystyle=$
$\displaystyle\sin^{2}(2\theta_{e\mu})\sin^{2}\left(\frac{\Delta
m^{2}_{41}L}{4E}\right)$ (15) $\displaystyle=$
$\displaystyle\sin^{2}(2\theta_{e\mu})\left(\frac{\Delta
m^{2}_{41}L}{4}\right)^{2}E^{-2}.$ (16)

The oscillation probability decreases as $E^{-2}$ for a fixed baseline. The
signal rates increase as $E^{3}$: there is a factor of $E^{2}$ from the solid
angle arising from the $1/\gamma$ opening angle and another factor of $E$ from
the cross section. The conclusion is that raising the stored muon energy will
increase the event rates linearly with energy for a fixed baseline. This
result has been confirmed by similar analyses for other muon-decay based
facilities (See sensitivity work in [4]).

##  9 Conclusion

The sensitivity of $\nu$STORM rules out the LSND 99% confidence interval at
$10\sigma$ using only appearance information. The appearance channel is the
CPT invariant of the observed anti-neutrino LSND anomaly. Optimizations have
been shown to guide future costing and performance work.

## Acknowledgements

The author thanks Alan Bross, John Cobb, and Joachim Kopp for their guidance
and knowledge. The author also thanks Ryan Bayes for the migration matrices
used in this analysis.

## References

  * [1] K.N. Abazajian, M.A. Acero, S.K. Agarwalla, A.A. Aguilar-Arevalo, C.H. Albright, et al.  Light Sterile Neutrinos: A White Paper.  2012\.  1204.5379.
  * [2] C. Athanassopoulos et al.  Evidence for nu(mu) —¿ nu(e) neutrino oscillations from LSND.  Phys.Rev.Lett., 81:1774–1777, 1998.  nucl-ex/9709006.
  * [3] Steve Baker and Robert D. Cousins.  Clarification of the use of chi-square and likelihood functions in fits to histograms.  Nuclear Instruments and Methods in Physics Research, 221(2):437 – 442, 1984.
  * [4] S. Choubey et al.  International Design Study for the Neutrino Factory, Interim Design Report.  2011\.  1112.2853.
  * [5] M. Aguilar-Benitez _et al._ Review of Particle Physics.  Physical Review D, 2008.
  * [6] Steve Geer, Olga Mena, and Silvia Pascoli.  A Low energy neutrino factory for large $\theta_{13}$.  Phys.Rev., D75:093001, 2007.  hep-ph/0701258.
  * [7] Carlo Giunti and Chung W. Kim.  Fundamentals of Neutrino Physics and Astrophysics.  Oxford University Press, USA, 2007.
  * [8] Carlo Giunti and Marco Laveder.  Towards 3+1 Neutrino Mixing.  2011\.  1109.4033.
  * [9] GNU.  GPL 3.0.
  * [10] M.C. Gonzalez-Garcia, Michele Maltoni, and Jordi Salvado.  Updated global fit to three neutrino mixing: status of the hints of theta13 ¿ 0.  JHEP, 1004:056, 2010.  1001.4524.
  * [11] Patrick Huber, Joachim Kopp, Manfred Lindner, Mark Rolinec, and Walter Winter.  New features in the simulation of neutrino oscillation experiments with GLoBES 3.0: General Long Baseline Experiment Simulator.  Comput.Phys.Commun., 177:432–438, 2007.  hep-ph/0701187.
  * [12] Patrick Huber, M. Lindner, and W. Winter.  Simulation of long-baseline neutrino oscillation experiments with GLoBES (General Long Baseline Experiment Simulator).  Comput.Phys.Commun., 167:195, 2005.  hep-ph/0407333.
  * [13] Joachim Kopp.  Efficient numerical diagonalization of hermitian $3\times 3$ matrices.  Int. J. Mod. Phys., C19:523–548, 2008.  physics/0610206.
  * [14] Joachim Kopp, Manfred Lindner, Toshihiko Ota, and Joe Sato.  Non-standard neutrino interactions in reactor and superbeam experiments.  Phys. Rev., D77:013007, 2008.  0708.0152.
  * [15] Mark D. Messier.  Evidence for neutrino mass from observations of atmospheric neutrinos with super-kamiokande.  1999\.  UMI-99-23965.
  * [16] $\nu$STORM Collaboration.  nuSTORM: Neutrinos from STORed Muons.  2012\.  1206.0294.
  * [17] E. A. Paschos and J. Y. Yu.  Neutrino interactions in oscillation experiments.  Phys. Rev., D65:033002, 2002.  hep-ph/0107261.
  * [18] Peter Renton.  Electroweak Interactions: An Introduction to the Physics of Quarks and Leptons.  Cambridge University Press, 1990.
  * [19] C. D. Tunnell.  https://code.launchpad.net/$\sim$c-tunnell1/+junk/, 5 2012.  questions should be directed to the corresponding author.
  * [20] Christopher D. Tunnell, John H. Cobb, and Alan D. Bross.  Sensitivity to eV-scale Neutrinos of Experiments at a Very Low Energy Neutrino Factory.  2011\.  1111.6550.
  * [21] Walter Winter.  Optimization of a Very Low Energy Neutrino Factory for the Disappearance Into Sterile Neutrinos.  2012\.  1204.2671.