Probing neutron star structure via f-mode oscillations and damping in dynamical spacetime models

TL;DR

This paper uses numerical relativity simulations to predict f-mode oscillation frequencies and damping times in neutron stars, providing a foundation for interpreting gravitational wave data to probe neutron star structure at high densities.

Contribution

It systematically studies the non-linear regime of neutron star f-modes using fully relativistic simulations, establishing accuracy benchmarks for future gravitational wave analyses.

Findings

F-mode frequencies match linear perturbation theory results.

Damping times deviate from linear predictions within numerical errors.

Provides blueprints for numerical accuracy in modeling neutron star oscillations.

Abstract

Gravitational wave and electromagnetic observations can provide new insights into the nature of matter at supra-nuclear densities inside neutron stars. Improvements in electromagnetic and gravitational wave sensing instruments continue to enhance the accuracy with which they can measure the masses, radii, and tidal deformability of neutron stars. These better measurements place tighter constraints on the equation of state of cold matter above nuclear density. In this article, we discuss a complementary approach to get insights into the structure of neutron stars by providing a model prediction for non-linear fundamental eigenmodes (f-modes) and their decay over time, which are thought to be induced by time-dependent tides in neutron star binaries. Building on pioneering studies that relate the properties of f-modes to the structure of neutron stars, we systematically study this link in the non-perturbative regime using models that utilize numerical relativity. Using a suite of fully relativistic numerical relativity simulations of oscillating TOV stars, we establish blueprints for the numerical accuracy needed to accurately compute the frequency and damping times of f-mode oscillations, which we expect to be a good guide for the requirements in the binary case. We show that the resulting f-mode frequencies match established results from linear perturbation theory, but the damping times within numerical errors depart from linear predictions. This work lays the foundation for upcoming studies aimed at a comparison of theoretical models of f-mode signatures in gravitational waves, and their uncertainties with actual gravitational wave data, searching for neutron star binaries on highly eccentric orbits, and probing neutron star structure at high densities.

Figures

Tables

Table 1. Table 1: This table lists all of the models in the first column and the important quantities that changed between models. These quantities were the spacetime formalism, the atmospheric density ρ atm subscript 𝜌 atm \rho_{\textrm{atm}} , the perturbation amplitude α 𝛼 \alpha , the Riemann solver, and the reconstruction method.

Model	Formalism	${\rho }_{atm}\left({\text{M}}^{-2}\right)$	$\alpha$	Riemann	Recon.
PPM Big $\alpha$	BSSN	${10}^{-10}$	0.141	Marquina	PPM
PPM Big ${\rho }_{\text{atm}}$	BSSN	${10}^{-10}$	0.0141	Marquina	PPM
PPM CCZ4	CCZ4	${10}^{-10}$	0.0141	Marquina	PPM
PPM	BSSN	${10}^{-14}$	0.0141	Marquina	PPM
WENO	BSSN	${10}^{-14}$	0.0141	Marquina	WENO
WENO HLLE	BSSN	${10}^{-14}$	0.0141	HLLE	WENO

Table 2. Table 2: Table depicting the results our from oscillating neutron star simulations. These results include the resolution of the simulation Δ ​ x Δ 𝑥 \Delta x , its frequency f 𝑓 f , the error in that frequency Δ ​ f Δ 𝑓 \Delta f , the damping time τ 𝜏 \tau and its estimated error Δ ​ τ Δ 𝜏 \Delta\tau . The damping error is the statistical fitting error of the data and does not include any systematic errors alluded to in Sections III.2 and III.6 . The table is divided by the simulation models described in Table 1 . Some cases yield convergence orders that exceed the formal accuracy of the numerical schemes involved. Such so-called overconvergence usually indicates that the data is close to but not quite within the convergent regime. At the end of each group of simulations, we include the damping time value at Δ ​ x → 0 → Δ 𝑥 0 \Delta x\rightarrow 0 using Richardson extrapolation. This procedure will systematically underestimate the error in presence of over-convergence and should be taken with a grain of salt. We did not include errors on this value nor the frequency when performing this extrapolation to infinitely fine resolution. The expected frequency and damping time for all simulations in the table according to linear theory are 1.579 kHz and 298 ms respectively.

	$\Delta x$ (M)	$f$ (kHz)	$\Delta f$ (kHz)	$\tau$ (ms)	$\Delta \tau$ (ms)	$n$
PPM Big $\alpha$	0.25	1.551	0.023	45.29	0.39	4.84
	0.16	1.552	0.013	58.02	0.93	4.84
	0.125	1.555	0.017	59.81	1.65	4.84
	0.08	1.555	0.013	63.76	1.42	4.84
	0.0625	1.555	0.011	64.39	1.46	4.84
	0	-	-	64.83	-	4.84
PPM Big ${\rho }_{\text{atm}}$	0.25	1.572	0.019	8.44	0.21	3.83
	0.16	1.579	0.014	38.36	1.56	3.83
	0.125	1.580	0.023	63.48	9.88	3.83
	0.08	1.580	0.004	178.73	3.79	3.83
	0.0625	1.581	0.004	262.11	7.58	3.83
	0	-	-	433.80	-	3.83
PPM CCZ4	0.25	1.574	0.004	8.41	0.13	3.14
	0.16	1.577	0.003	35.36	0.96	3.14
	0.125	1.579	0.004	73.21	0.57	3.14
	0.08	1.580	0.004	195.57	3.76	3.14
	0.0625	1.581	0.004	268.25	6.33	3.14
	0	-	-	376.65	-	3.14
PPM	0.25	1.573	0.005	7.47	0.16	2.98
	0.16	1.576	0.004	36.61	0.55	2.98
	0.125	1.576	0.003	73.35	1.05	2.98
	0.08	1.577	0.004	168.48	2.71	2.98
	0.0625	1.577	0.003	222.20	3.85	2.98
	0	-	-	313.74	-	2.98
WENO	0.25	1.573	0.005	80.61	0.99	2.79
	0.16	1.576	0.004	165.98	1.54	2.79
	0.125	1.577	0.004	209.60	2.14	2.79
	0.08	1.577	0.003	258.09	4.07	2.79
	0.0625	1.578	0.003	270.44	4.80	2.79
	0	-	-	283.80	-	2.79
WENO HLLE	0.25	1.574	0.005	124.48	8.42	5.77
	0.16	1.576	0.004	195.72	1.71	5.77
	0.125	1.577	0.003	232.12	2.44	5.77
	0.08	1.579	0.011	265.56	2.80	5.77
	0.0625	1.578	0.003	275.47	5.00	5.77
	0	-	-	290.37	-	5.77

Equations

$$R_{\mu \nu }-\frac{1}{2}R{g}_{\mu \nu }=8\pi T_{\mu \nu }$$

$$d{s}^2=-{\alpha }^{2}d{t}^2+{\gamma }_{ij}(d{x}^i+{\beta }^{i}dt)(d{x}^j+{\beta }^{j}dt)$$

$$n^{\mu }=\left({\alpha }^{-1},-\frac{{\beta }^i}{\alpha }\right)$$

$$8\pi \left(T_{\mu \nu }-\frac{1}{2}{g}_{\mu \nu }T\right)$$

$$Z_{\mu }=0$$

$$T_{\mu \nu }=\rho h{u}_{\mu }{u}_{\nu }+p{g}_{\mu \nu }$$

$${\nabla }_{\mu }{T}^{\mu \nu }=0,$$

$${\nabla }_{\mu }\left(\rho u^{\mu }\right)=0$$

$${\partial }_t\mathbf{U}+{\partial }_i{\mathbf{F}}^i=\mathbf{S}$$

$$\frac{\boldsymbol{U}}{\sqrt{\gamma}}=\left(\begin{array}[]{c}D\\ S_{j}\\ \tau\end{array}\right)=\left(\begin{array}[]{c}\rho W\\ \rho hW^{2}v_{j}\\ \rho hW^{2}-p-\rho W\end{array}\right)$$

$$\frac{\boldsymbol{F}^{i}}{\sqrt{\gamma}}=\left(\begin{array}[]{c}D\left(\alpha v^{i}-\beta^{i}\right)\\ \alpha\tilde{S}_{\,j}^{i}-\beta S_{j}\\ \alpha\left(S^{i}-Dv^{i}\right)-\tau\beta^{i}\end{array}\right)$$

$$\delta P=\alpha \rho {\left(\frac{r}{R}\right)}^2{Y}_{22}$$

$$P$$

$${\ddot{h}}_{+}-\mathrm{i}{\ddot{h}}_{\times }={\psi }_4=\sum _{\ell =2}^{\infty }\sum _{m=-\ell }^{\ell }{\psi }_4^{\ell m}{}_{-2}{Y}_{\ell m}(\theta ,\phi ),$$

$${\psi }_{4,\text{high}}-{\psi }_{4,\text{med}}=c_n({\psi }_{4,\text{med}}-{\psi }_{4,\text{low}}),$$

$$c_n=\frac{\Delta x_{\text{high}}^n-\Delta x_{\text{med}}^n}{\Delta x_{\text{med}}^n-\Delta x_{\text{low}}^n}.$$

$$\text{error}=|c_n\Delta x^n|.$$

Full text

I Introduction

Neutron stars Landau (1932) present unique opportunities to study the nature of matter at supra-nuclear densities at relatively low temperatures, and in the presence of extreme gravitational interactions. Given the diverse scenarios that have been invoked to describe the properties of cold, dense matter in the cores of neutron stars (nucleons, hyperons or free quarks) Stone (2016); Alcock et al. (1986); Yakovlev and Pethick (2004); Weber (2005), it is essential to confront these theories with observations in the gravitational or electromagnetic spectrum, or a combination of them through multimessenger astrophysics discovery campaigns Abbott et al. (2017a, b). These studies will constrain the space of permissible models and hopefully, gain a better understanding of the physics at supranuclear densities. Electromagnetic observations of the pulsars PSR J1614-2230 Fonseca et al. (2016) and PSR 0348+0432Antoniadis et al. (2013), which are consistent with neutron stars with masses $M\sim 2M_{\odot}$, have already significantly reduced the existing number of astrophysically viable equations of state (EOS), $P(\epsilon)$, i.e., the relation between the pressure $P$ and the total energy density $\epsilon$. In Rezzolla et al. (2018) an upper limit for the neutron star maximum mass of $M=2.16M_{\odot}$ was derived based on GW170817 and basic arguments on kilonovas. A consistent limit was derived from different considerations Margalit and Metzger (2017) base on previous arguments about the nature of short gamma-ray bursts Bauswein et al. (2013); Fryer et al. (2015); Lawrence et al. (2015).

X-ray observations of the pulse profiles generated by hot spots on moderately spinning neutron stars with the ongoing NICER Gendreau et al. (2012) and future LOFT Feroci et al. (2012) missions may enable measurements of the neutron star radii at the 5% level Psaltis et al. (2014), improving existing measurements that are dominated by systematic errors Miller (2013); Miller and Lamb (2016). This precision can be achieved through long exposure times using NICER, or using the large collecting area of LOFT.

LIGO Abbott et al. (2016); The LIGO Scientific Collaboration et al. (2015) and Virgo Acernese et al. (2015) gravitational wave (GW) observations of the neutron star merger GW170817 Abbott et al. (2017a, b); Coulter et al. (2017) imply that the two neutron stars had radii $R_{\{1,\,2\}}=\tensor{*}[]{\mathrm{11.9}}{{}^{+1.4}_{-1.4}}\,\textrm{km}$, and to place constraints on the EOS at supra-nuclear densities $P(\epsilon)=\tensor{*}[]{\mathrm{3.5}}{{}^{+2.7}_{-1.7}}\times 10^{34}\,\textrm{dyn}\,\textrm{cm}^{-2}$ (quantified at the 90% confidence level The LIGO Scientific Collaboration et al. (2018)) as well as the dimensionless tidal deformability \tilde{\Lambda}\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower 4.0pt\hbox{\sim}}\hss}\hbox{<}}}800 Abbott et al. (2017a). As a reference point, the GW measurement of $P(\epsilon)$ for GW170817 is at twice nuclear saturation density, whereas terrestrial laboratories can only test and constrain cold EOSs at densities near or below saturation density of nuclei. GW observations enable these measurements because the waveform models used to analyze LIGO and Virgo data encode information about the static and dynamic tides that determine the dynamical evolution of the system, driving it to merger faster than a system of two point masses Blanchet (2006); Rosswog et al. (2003); Haas et al. (2016); Parisi and Sturani (2018).

GW and electromagnetic observations will continue to provide new and detailed information about the nature of neutron star matter. GW observations will shed light on the masses and tidal deformability of neutron stars, while electromagnetic observations will enable accurate measurements of their mass and radius Gendreau et al. (2012) and may even be informative for studies of stellar structure in the context of non-GR spacetimes Mendes and Ortiz (2018).

Even though multimessenger observations will provide new and detailed information about the properties of neutron stars, they are likely to leave open questions, especially in light of potential emission-model dependence, which calls for independent and complementary measurements. In this article, we will focus on a different observational scenario to obtain this information from direct observations of the stars, namely, the measurement of fundamental eigenmodes (f-modes) of the individual stars, which can be caused by tidal effects, and for which pressure is the restoring force.

Proposed scenarios for the generation of f-modes in neutron stars include neutron star binaries or neutron star-black hole systems in highly eccentric orbits Gold et al. (2012); East et al. (2012); East and Pretorius (2012); East et al. (2016a, b); Papenfort et al. (2018), i.e. binary separations comparable to the radii of the binary components, such that violent changes in the tidal field can perturb the stars without necessarily causing an immediate collision Chirenti et al. (2017); Huerta et al. (2017, 2014); Huerta and Brown (2013); Huerta et al. (2018); Rebei et al. (2018). These close interactions will induce f-mode oscillations in cold, slowly rotating neutron stars, thereby simplifying the analysis and modeling of f-modes as compared to the differentially rotating scenario Chirenti et al. (2017) or shock heated matter in hypermassive neutron stars (HMNSs) Bauswein et al. (2014); Takami et al. (2014).

Extracting physics encoded in the f-modes from observational data involves an accurate model prediction Steinhoff et al. (2016); Yang et al. (2018), and measurement of their frequencies $\omega$ and damping times $\tau$. Seminal studies have established that these two observables are related to the average density and compactness Andersson and Kokkotas (1998); Benhar et al. (2004). In Tsui and Leung (2005); Lau et al. (2010) universal relations using the effective compactness $\eta=\left(M/I^{3}\right)^{1/2}$ with mass $M$ and moment of inertia $I$ were discovered, and were tested in Chirenti et al. (2015).

In principle, once $\omega$ and $\tau$ are extracted from GW observations, it may be possible to infer $M$ and $I$. Furthermore, the amplitude of the f-mode depends on the orbital energy deposited in the mode, and this can be used to measure the Love number $\lambda$ and the rotationally-induced quadrupole $Q$ through $I$-Love-$Q$ relations Flanagan and Hinderer (2008); Yagi and Yunes (2013a, b, 2014).

In this article, we carry out a detailed analyses to quantify the required accuracy in numerical relativity simulations to extract the f-mode oscillations and damping times. To do this, we simulate oscillating neutron stars using the open source Einstein Toolkit numerical relativity software, and compare these results to linear perturbation theory. This work lays the foundation for future numerical relativity simulations of highly eccentric neutron star mergers (see Dietrich et al. (2018) for a state-of-the-art library) with the aim of accurately extracting f-mode oscillation damping times and frequencies as well as a systematic investigation for a variety of equations of state and the effect of model uncertainties on parameter estimation. A successful detection and accurate parameter estimation from eccentric neutron star binaries would provide a wealth of information on astrophysical evolution scenarios Gondán and Kocsis (2018).

This article is organized as follows: Section II outlines the approach we have followed to couple Einstein’s field equations with a relativistic perfect fluid. We describe the construction of initial data to excite f-mode oscillations in Tolman-Oppenheimer-Volkoff (TOV) stars, and the extraction of the corresponding GWs. In Section III, we show that the numerical formalism introduced in this study exhibits the necessary convergence to accurately compute the frequency and damping time of f-mode oscillations. We also demonstrate that our simulations reproduce the expected results from linear perturbation theory in the appropriate limit. We summarize our findings and future directions of work in Section IV.

II Methods

We evolve the spacetime by solving the Einstein field equations using the tools of numerical relativity.

[TABLE]

We express the four-dimensional metric $g_{\mu\nu}$ in the standard 3+1 split of spacetime Alcubierre (2008); Bona et al. (2009); Baumgarte and Shapiro (2010); Gourgoulhon (2012); Rezzolla and Zanotti (2013); Shibata (2016) as

[TABLE]

with the spatial metric $\gamma_{ij}$ induced on each hypersurface of the spacetime foliation, the lapse $\alpha$, shift $\beta^{i}$, and the hypersurface unit normal vector

[TABLE]

Projecting the Einstein equations onto the hypersurface and along $n^{\mu}$ yields a set of constraint and evolution equations that together with suitable initial data form our initial value problem.

II.1 Spacetime Evolution

There are infinitely many formulations of the Einstein equations that are equivalent in the continuum limit of infinite spatial and temporal resolution, but that have distinct principal parts and hyperbolicity properties, and hence behave differently at finite resolution. In this study we primarily employed the BSSN formalism Shibata and Nakamura (1995); Baumgarte and Shapiro (1999).

The CCZ4 formalism Alic et al. (2012) was also used for some results. The starting point is the covariant extension of the Einstein equations

[TABLE]

also known as the Z4 formulation, with the four-vector $Z_{\mu}$ containing the four constraints. These constraints have to be satisfied for all times and in particular by the initial data used as input for the time integration. If

[TABLE]

then the Einstein field equations are recovered. This equivalence in the continuum limit implies that solutions to Eq. (II.1) under the constraint (5) are also solutions to Eq. (1). The additional terms involving the constraint damping coefficients $\kappa_{1,2}$ lead to an exponential damping of the constraint violations on a characteristic timescale proportional to the coefficients. Even in the zero damping case $\kappa_{1}=\kappa_{2}=0$, Z4-based formulations feature a propagating mode associated with the Hamiltonian constraint violations as opposed to the classic BSSN formulation Shibata and Nakamura (1995); Baumgarte and Shapiro (1999), where a 0-speed mode leads to local build-up of Hamiltonian constraint violations due to numerical error. In Z4-based formulations such violations propagate and therefore have a chance to leave the grid. It is therefore expected that Z4-based formulations yield more accurate spacetime evolutions than the BSSN system. Employing a conformal and traceless decomposition of the resulting equations leads to the CCZ4 system, see Alic et al. (2012) for details and Hilditch et al. (2013); Bernuzzi and Hilditch (2010) for an alternative approach. For the cases in this work in which the CCZ4 formulation was used, we set $\kappa_{1}=0.02$, $\kappa_{2}=0$.

With that said, the CCZ4 formulation is sensitive to the boundary conditions and may lead to instabilities if not applied properly Alic et al. (2012). To ensure robustness for our work, BSSN was employed for most of this work.

The evolution of the gauge conditions, i.e., lapse and shift, is governed by standard 1+log slicing condition and a Gamma driver (see e.g.Baumgarte and Shapiro (2010); Baiotti and Rezzolla (2017)).

We are making use of the publicly available McLachlan code Brown et al. (2009) which implements the above evolution equations with fourth order central finite differences and is part of the Einstein Toolkit Löffler et al. (2012).

The non-linear stability of smoothness properties of the spacetime evolution is enhanced by adding an artificial Kreiss-Oliger dissipation Kreiss and Oliger (1973) to the spacetime variables.

II.2 General Relativistic Hydrodynamics

The neutron stars are modeled by a relativistic perfect fluid

[TABLE]

where $u^{\mu}$ is the 4-velocity, $\rho$ is the rest mass, $h$ is the specific enthalpy and $p$ is the pressure of the fluid.

The standard relativistic hydrodynamic equations are the local conservation law for the energy-momentum tensor

[TABLE]

and the conservation of rest mass

[TABLE]

Both equations together with the spacetime metric $g_{\mu\nu}$ represent the equations of motion of the fluid, which are closed by a suitable equation of state (EOS) with the general form $p=p\left(\rho,T,Y_{e}\right)$ where $T$ is the temperature and $Y_{e}$ the so-called electron fraction (low $Y_{e}$ indicates neutron-rich material). In this pilot study we restrict ourselfs to polytropes, but in the future we intend to investigate more realistic equations of state.

These equations can be cast into conservative form, which is also known as the “Valencia formulation” Banyuls et al. (1997)

[TABLE]

where the conserved variables $\boldsymbol{U}$ are defined by

[TABLE]

and the fluxes as well as the sources are given by

[TABLE]

These conservation equations are implemented in the Einstein Toolkit’s GRHydro code Mösta et al. (2014), which uses high-resolution shock-capturing methods.

We employ a fourth order Runge-Kutta method to advance in time, either the WENO or PPM reconstruction method, second order flux evaluations and either the HLLE or Marquina Riemann solver. We set the Courant factor to 0.4.

Since we are not able to accurately track regions that transition from the fluid to the vacuum regime, all simulations have an artificial low-density background atmosphere which is evolved freely. For all simulations presented in this study we choose an atmosphere value of $\rho_{\textrm{atm}}=10^{-10}\text{M}^{-2}\approx 6.2\times 10^{7}\textrm{g/cm}^{3}$ or $\rho_{\textrm{atm}}=10^{-14}\,\text{M}^{-2}\approx 6.2\times 10^{3}\textrm{g/cm}^{3}$.

II.3 Initial Data

We use a perturbed TOV solution as initial data. The matter perturbation is implemented by adding the (2,2) pressure eigenfunction of an incompressible Newtonian star in the Cowling approximation Jones et al. (2002). This form of the perturbation is smooth and provides a well controlled framework to quantify numerical accuracy.

Specifically, the initial data was created with the Einstein Toolkit’s TOVSolver thorn. Then, we added a pressure perturbation to excite this f-mode oscillations using the Cowling approximation. The perturbation was given by

[TABLE]

where $\delta P$ is the perturbed pressure, $\alpha$ is the amplitude of the perturbation, $\rho$ is density, $r$ is radial distance from the center of the star, $R$ is the radius of the star, and $Y_{22}$ is the $l=2$ $m=2$ spherical harmonic. A plot of the density distribution resulting from this perturbation can be seen in Fig 1.

II.4 Equation of state

We adopt a polytropic equation of state

[TABLE]

where $\Gamma$ is the adiabatic index and $K$ is the constant of proportionality. For this work, we examined a polytrope with $\Gamma=2$ and $K=100$ which is a common choice in the literature. We leave to future work the expansion of the methods applied in this paper to other EOSs.

II.5 Grid Setup

The numerical grid is managed by the mesh refinement Carpet Schnetter et al. (2004) driver for Cactus Goodale et al. (2003).

It implements a non-uniform grid via a nested set of movable boxes (box-in-box) together with a hierarchical (Berger-Oliger-style) timestepping. In these simulations we employ four refinement levels where each additional level doubles the resolution of the enclosing one. On the finest grid, which is separately centered around each neutron star, the resolution in each dimension is $\Delta x=0.0625\text{M}\approx 93.75\rm{m}$ for the highest resolution. The physical domain extends to $128\text{M}\approx 192\rm{km}$.

II.6 Gravitational Waves

The gravitational wave strain is extracted using the standard Newman-Penrose formalism Newman and Penrose (1962), in which a particular contraction of the Weyl tensor with a suitably chosen null tetrad yields a gauge-invariant quantity $\psi_{4}$ that encodes the outgoing gravitational radiation Bishop and Rezzolla (2016). Specifically, $\psi_{4}$ is related to the second time derivative of the two strain polarizations $\ddot{h}_{+,\times}$ by

[TABLE]

where we also introduced the multipole expansion of $\psi_{4}$ in spin-weighted spherical harmonics Goldberg et al. (1967) of spin-weight $s=-2$.

II.7 Models

We employed various models in which we altered simulation parameters such as perturbation amplitude, the spacetime formalism, etc. We present and label these models in Table 1. All models used a $1.4M_{\odot}$ TOV star which had a central density of $0.00128\text{M}^{-2}$ and the EOS described in Section II.4.

III Results

III.1 Extracted Gravitational Waveforms

We extracted the gravitational waves at distances of 10M, 40M, 70M, 100M, and 120M. The analysis relied on the $l=2$ $m=2$ mode of the $\psi_{4}$ data extracted from 100M unless otherwise stated.

This system produced linearly polarized gravitational waves, which differs from other phenomena in gravitational wave physics which produce elliptically polarized waves. This linear polarization arises from the lack of rotation in our oscillating neutron star. In contrast, more often studied gravitational wave producing systems such as binary neutron star mergers and binary black hole mergers, use their rotation to produce their gravitational waves and give their gravitational waves their elliptical polarization.

The polarization is imprinted in the real and imaginary parts of $\psi_{4}$ which represent the plus and cross polarizations respectively. The real part exhibits a sinusoidal form with a frequency close to the fundamental oscillation mode frequency and on top a secular exponential decay, see Fig 2. The imaginary part is zero to almost machine precision, indicating that no cross polarization was present.

We note that $\psi_{4}$ had high frequency noise as a result from passing through refinement layers. We removed these high frequency features with low pass filters that are described in detail in Appendix A.

III.2 Convergence

We demonstrate that differences between resolutions behave as expected given the numerical method we adopt. This implies that the numerical error in our simulations can be quantified. Only with such a reliable error estimate is it possible to tell whether any damping in the gravitational waves is physically caused by gravitational radiation or by numerical dissipation. There are several possibilities to demonstrate convergence in the GW output of such models. We choose to perform convergence tests on the extrema of $\psi_{4}$.

As usual, we assume that the numerical error as represented by a Taylor expansion is dominated by the leading order term. In this case three resolutions suffice to demonstrate convergence and fit the parameters of the Taylor expansion. Specifically, the numerical solution for high, medium, and low resolutions are related by

[TABLE]

where $c_{n}$ is the convergence factor of order $n$ and is determined only by the resolutions and convergence order as

[TABLE]

We present the convergence order for all our numerical relativity simulations in Table 2. A sample case of this analysis for the WENO reconstruction is presented in Fig 3. Using the convergence factor and order, we were able to estimate the systematic error due to the changing resolution as

[TABLE]

The results for the estimated error of the WENO reconstruction is presented in Fig 4. This error behaves as systematic error caused by the differences in simulation resolution. As such, we did not propagate this error into the analysis because the analysis aims at characterizing the effect of this same resolution dependence on the damping times.

III.3 Fourier Analysis

Frequency spectra of the gravitational waveform were examined by taking the Fast Fourier Transform (FFT) of the $\psi_{4}$ data. We plot the frequency spectrum of one of the runs in Fig 5. We see that the most excited frequency is nearly aligned with the expected frequency of 1.579 kHz, illustrating that our perturbation excited the desired f-mode. In addition, we observe at least three secondary excitations which we believe to be p-modes.

III.4 Frequency

We also investigated the frequency evolution as the run progressed. Therefore, we took the roots of $\psi_{4}$ time series to estimate the period, which was in turn used to obtain frequency estimates for each cycle. We used the average frequency as the value for the simulation and the standard deviation as the error.

Although most of the runs exhibited constant frequency regardless of resolution, runs with high density atmospheres $10^{-10}\text{M}^{-2}$ of and a length longer than 10 ms experienced a rise in frequency as runs progressed.

In contrast, the frequency of low density atmosphere runs of $10^{-14}\text{M}^{-2}$ remained constant for even the longest runs which laster 50 ms. We can see the frequency evolution in Fig 6. The likely physical cause of this behavior is accretion of artificial atmosphere material onto the star which leads to a steady change in the background stellar model and hence a slight change in its fundamental oscillation frequencies.

We found that the frequency was fairly independent of resolution as displayed in Table 2. Though the higher amplitude runs did experience a lower frequency than the lower perturbation runs. The frequencies all appear to be around 1.579 kHz which is the expected frequency for f-mode oscillations based on linear theory.

We note that since the gravitational waves were linearly polarized, we could not determine instantaneous frequency by taking the derivative of the phase of the gravitational waves as the imaginary part was nearly zero.

III.5 Damping Times

Damping times were extracted by performing linear fits to the natural log of the extrema of $\psi_{4}$, as shown in Fig 7. Error estimates were computed using the standard deviations of the damping time calculations using the extrema of $\psi_{4}$ individually rather than from all the extrema in the simulation dataset.

To compute the damping times from our suite of numerical relativity simulations we used the WENO and PPM reconstructions methods, finding that the former is much more robust and accurate than the latter. Since simulations with $\Delta x<0.0625\text{M}$ were deemed to be too computationally intensive to run, we extrapolated these results to the continuum limit to see if they approach the expected value of 298 ms. We used Richardson extrapolation with unequal resolution ratios to estimate the value at infinitely fine resolution. A summary of these calculations is presented in Table 2.

Fig 8 shows that the WENO reconstruction approaches the linear perturbation prediction of 298 ms for the highest resolution run of each simulation. The three simulations that significantly depart from the linear perturbation regime correspond to two categories: (i) the two simulations with damping times much larger than 298 ms correspond to simulations with large artificial atmospheres, $\rho_{\textrm{atm}}\sim 10^{-10}\text{M}^{-2}$; (ii) the simulation with estimated damping time much smaller than 298 ms is clearly outside of the linear perturbation regime by construction.

These results indicate that provided a numerical setup in which the atmosphere used for the simulations satisfies \rho_{\textrm{atm}}\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower 4.0pt\hbox{\sim}}\hss}\hbox{<}}}10^{-14}\text{M}^{-2}, and the perturbation applied to the TOV star is within the linear regime, the extracted damping times with the WENO reconstruction method from our suite of numerical simulations is consistent with linear perturbation predictions.

III.6 Evolution of Damping Time

High resolution runs with longer damping times required longer analysis times to obtain satisfactory statistical fits. This necessity derives from the difference in position of individual points being significant compared to the difference induced by the true damping of the waveform over short timescales. In essence, the longer we run the simulation, the closer we get to running for the damping time we would like to measure; thereby, we reduce the statistical fitting error in the damping time.

Although increasing the run time decreases the statistical error in the damping time, it introduces systematic error in the damping time calculation. We found that as the run progressed, the waveform did not exhibit a true exponential decay. Instead, the waveform systematically decayed quicker than an exponential damping time would suggest at long times. This posits the existence some other source of decay in the waveforms, most likely some form of numerical dissipation.

Fig 9 illustrates the evolution of the damping time value vs the analysis time, the cutoff time for waveform data to in the analysis to calculate the damping time beyond which the data would not be used for such analysis. We also note that the first 2 ms of data also was not used in the analysis to allow time for any transients to die out. In addition, the figure depicts the statistical fitting error at each of those points. One can observe that the damping time and statistical error decrease as the analysis time increase. We observed a local minimum in the statistical error at an analysis time of 15 ms and thus chose that analysis time for the computation of the damping times in Table 2.

IV Conclusions

We have carried out systematic analyses of the accuracy needed to extract the frequency and damping time of f-mode oscillations. We have demonstrated the robustness of our results by extracting these observables, and demonstrating that our results are similar when using both the WENO and PPM reconstruction methods.

We have shown that for small perturbations, the damping times extracted from our numerical relativity simulations are consisted with linear theory perturbations when we consider atmospheres $\rho_{\textrm{atm}}\leq 10^{-14}\text{M}^{-2}$. On the other hand, we have found that if artificial atmospheres are poorly chosen, $\rho_{\textrm{atm}}\leq 10^{-10}\text{M}^{-2}$, the extracted damping times from our numerical simulations differ significantly from linear perturbation predictions at evolutions \mathrel{\hbox{\hbox to0.0pt{\hbox{\lower 4.0pt\hbox{\sim}}\hss}\hbox{>}}}10ms. This may be caused due to accretion of the atmosphere, which introduces a secular change (increased mass) in the background solution.

The convergence analyses we present suggest that for a careful choice of numerical schemes a minimum resolution of $\sim 0.125\text{M}$ is needed to accurately extract the GW induced damping time of f-mode oscillations. WENO reconstruction clearly outperforms the PPM scheme, which needed twice the resolution to obtain similar accuracies in the damping times. This resolution is within the range of existing neutron star numerical relativity simulations that have been presented in the literature albeit on the more computationally expensive end.

Having established a blueprint of numerical accuracy to extract the frequency and damping times of f-modes, in future work we will put this formalism to use and investigate the accuracy with which we can infer properties of the EOS through GW observations of highly eccentric neutron star encounters.

Acknowledgements

This research is part of the Blue Waters sustained-petascale computing project, which is supported by the National Science Foundation (awards OCI-0725070 and ACI-1238993) and the State of Illinois. Blue Waters is a joint effort of the University of Illinois at Urbana-Champaign and its National Center for Supercomputing Applications. Some numerical results were obtained using the University of Maryland Deepthought2 HPC Cluster. This work was supported in part by joint research workshop award 2015/50421-8 from FAPESP and the University of Maryland, College Park. This work was supported in part by Perimeter Institute for Theoretical Physics. Research at Perimeter Institute is supported by the Government of Canada through the Department of Innovation, Science and Economic Development Canada and by the Province of Ontario through the Ministry of Research, Innovation and Science. We acknowledge support from the NCSA. NSF-1550514 and NSF-1659702 grants are gratefully acknowledged. We thank the NCSA Gravity Group for useful feedback.

Appendix A Filtering

Initially, the $\psi_{4}$ data contained significant high frequency noise that interfered with our analysis. We noted that the signals extracted from 10M, before passing through the first refinement boundary, did not exhibit any such high frequency noise. Therefore,we hypothesized that passing through the various refinement boundaries as it traveled to the relevant extraction sphere at 100M.

To compensate for this high frequency behavior induced from the refinement layers, we applied a low pass filter to the raw data. Specifically, we applied a butter filter with zero phase filtering. Relative to the sample frequency, the filter had a passband frequency of 0.01 and a stopband frequency of 0.04. The passband ripple amplitude of the filter was 1 dB and the stopband attenuation was 60 dB. This filtering was essential to ensuring the extracted $\psi_{4}$ produced meaningful data analysis especially for the damping times. Figs 10 and 11 illustrates the effect of the filtering on the data.

We observe the frequency spectrum before and after the filtering by looking at the looking at the FFTs extracted at 100M in Fig 12. We observe that the unfiltered frequency spectrum is very noisy, which is most apparent after the second peak. The filtering eliminates this noise. However, the filtering also significantly suppresses the second peak and renders the third and fourth invisible.