# Axially symmetric pseudo-Newtonian hydrodynamics code

Jinho Kim,1 Hee Il Kim,1,211footnotemark: 1 Matthew William
Choptuik,3,411footnotemark: 1 and Hyung Mok Lee111footnotemark: 1
1Department of Physics and Astronomy, FPRD, Seoul National University, Seoul,
151-742, Korea
2Korea Institute of Science and Technology Information, 245 Daehak-ro,
Yuseong-gu, Daejeon, 305-806, Korea
3Department of Physics and Astronomy, University of British Columbia,
Vancouver, Canada
4Canadian Institute for Advanced Research Cosmology & Gravity Program E-mail:
<EMAIL_ADDRESS>(JK<EMAIL_ADDRESS>(HIK<EMAIL_ADDRESS>(MWC<EMAIL_ADDRESS>(HML)

(Released 2012 Xxxxx XX)

###### Abstract

We develop a numerical hydrodynamics code using a pseudo-Newtonian formulation
that uses the weak field approximation for the geometry, and a generalized
source term for the Poisson equation that takes into account relativistic
effects. The code was designed to treat moderately relativistic systems such
as rapidly rotating neutron stars. The hydrodynamic equations are solved using
a finite volume method with High Resolution Shock Capturing (HRSC) techniques.
We implement several different slope limiters for second order reconstruction
schemes and also investigate higher order reconstructions such as PPM, ENO and
WENO. We use the method of lines (MoL) to convert the mixed spatial-time
partial differential equations into ordinary differential equations (ODEs)
that depend only on time. These ODEs are solved using second and third order
Runge-Kutta methods. The Poisson equation for the gravitational potential is
solved with a multigrid method, and to simplify the boundary condition, we use
compactified coordinates which map spatial infinity to a finite computational
coordinate using a tangent function. In order to confirm the validity of our
code, we carry out four different tests including one and two dimensional
shock tube tests, stationary star tests of both non-rotating and rotating
models and radial oscillation mode tests for spherical stars. In the shock
tube tests, the code shows good agreement with analytic solutions which
include shocks, rarefaction waves and contact discontinuities. The code is
found to be stable and accurate: for example, when solving a stationary
stellar model the fractional changes in the maximum density, total mass, and
total angular momentum per dynamical time are found to be $3\times 10^{-6}$,
$5\times 10^{-7}$ and $2\times 10^{-6}$, respectively. We also find that the
frequencies of the radial modes obtained by the numerical simulation of the
steady state star agree very well with those obtained by linear analysis.

###### keywords:

relativistic processes - gravitation - hydrodynamics hydrodynamics- methods:
numerical

††pagerange: Axially symmetric pseudo-Newtonian hydrodynamics code–A††pubyear:
2012

##  1 Introduction

It is necessary to take into account both special- and general relativistic
effects in the studies of the dynamics of compact astrophysical object such as
neutron stars and black holes. Some pulsars produce pulses of up to 1 KHz,
corresponding to rotation speeds at the surface of around $0.2c$. Their
typical sizes and masses are known to be around $10\rm{km}$ and $1.4\sim
2\mbox{$\rm{M}_{\odot}$}$, respectively, giving compactness,
$GM/Rc^{2}=0.2\sim 0.3$. Therefore, a Newtonian approach cannot properly
describe neutron stars, even for the non-rotating case.

In general relativity, the dynamics of gravity (or spacetime) can be studied
by solving the Einstein equations. The equations of motion for the matter are
given, in part, by the conservation law of the energy-momentum tensor which
itself sources the gravitational field. Computational approaches for solving
general relativistic field equations constitute the field of numerical
relativity.

Over the past few decades, many general relativistic hydrodynamic codes have
been developed, starting with Wilson (1972) who proposed a 3+1 Eulerian
formulation (see also Wilson & Mathew, 2003). Although Wilson’s numerical
approach was widely used to study problems such as core collapse and accretion
disks, it produced large errors when fluid flows became ultra-relativistic
(Centrella & Wilson, 1984; Norman & Winkler, 1986). In order to avoid these
excessive errors, a new formulation was proposed by Marti et al. (1991). This
formulation makes it possible to use existing numerical techniques based on
characteristic approaches for Newtonian hydrodynamics. In particular, these
include High Resolution Shock Capturing (hereafter HRSC) methods that reduce
the order of accuracy near shocks, but minimize the amount of numerical
dissipation. This dissipation is very unnatural and can result in non-physical
effects in the numerical results. Marti’s formulation was extended to the
general relativistic case by the Valencia group (Font et al., 2000), and this
last work forms the basis for most recent general relativistic hydrodynamical
codes. Recent reviews of the formulation and numerical methods can be found in
Martí & Müller (2003) and Font (2008).

However, when working in multiple spatial-dimensions, it still requires a lot
of computational resources to treat fluid dynamics in concert with the
evolution of the general relativistic gravitational field. In addition,
numerical relativity simulations have frequently encountered instabilities
which are often associated with violations of the Hamiltonian and momentum
constraints. (However, with the development of new formulations which cast the
Einstein equations in appropriate hyperbolic forms, as well as the use of
constraint-damping techniques, significant progress has been made on this
front: see Sarbach & Tiglio, 2012 for a very recent review of this subject).
For these reasons, simulations using Newtonian gravity are still used even
though they are not applicable to very compact objects.

The aims of this paper are 1) to introduce a new formulation which applies a
pseudo-Newtonian approach (Kim et al., 2009) to the study of moderately
relativistic objects and 2) to describe a numerical implementation of this
method. In our pseudo-Newtonian approach, which was introduced by Kim et al.
(2009) for steady state models, the gravitational field is treated by a weak
field approximation, but special relativistic effects are correctly taken into
account. Specifically, the Newtonian gravitational potential that appears in
the weak field metric satisfies a Poisson equation, but the mass density that
appears as a source term for that equation is modified to include relativistic
effects. Of course this method cannot be applied to highly relativistic
systems, but Kim et al. (2009) showed that the pseudo-Newtonian formulation is
valid for the modeling of mildly compact objects, such as rotating neutron
stars having surface rotation velocity up to $\sim 0.2c$ and compactness $\sim
0.2$ (Kim et al., 2009). In this paper, we extend the pseudo-Newtonian
approach to hydrodynamical systems where the flows can be ultra-relativistic
and gravity can be moderately strong.

The remainder of the paper is structured as follows: In section 2, we present
the formulation and governing equations for our system, while the numerical
techniques employed in our study are given in section 3. We discuss various
numerical tests of our code’s treatment of hydrodynamics for the case of shock
tubes in 4, and for stationary stars in 5. A test which compares radial
pulsation mode frequencies for polytropic stars determined through dynamical
evolution to those computed in linear theory is detailed in section 6. We
conclude with a summary and discussion in section 7.

Throughout this paper we use units in which $c=G=\mbox{$\rm{M}_{\odot}$}=1$:
these correspond a to unit time $=4.92\times 10^{-3}\rm{ms}$, unit length
$=1.47\,\rm{km}$ and unit mass $=1.99\times 10^{33}\rm{g}$.

##  2 Formulation

Our pseudo-Newtonian method was first discussed in the steady state context by
Kim et al. (2009). We assume the weak field metric,

$\begin{split}ds^{2}&=g_{\mu\nu}dx^{\mu}dx^{\nu}\\\
&=-(1+2\Phi)dt^{2}+(1+2\Phi)^{-1}\delta_{ij}dx^{i}dx^{j},\end{split}$ (1)

where $g_{\mu\nu}$ is the spacetime metric and $\Phi$ is the Newtonian
gravitational potential. With this metric, we neglect all higher order effects
such as frame dragging and describe gravity using only a single gravitational
potential, just as in the Newtonian case. The gravitational potential
satisfies a Poisson equation with the active mass density,
$\rho_{\rm{active}}$, providing the source:

$\nabla^{2}\Phi=4\pi\rho_{\rm{active}}.$ (2)

The active mass density is computed from the relativistic definition of the
energy. For a perfect fluid, the energy momentum tensor can be expressed as

$T^{\mu\nu}=\rho_{0}hu^{\mu}u^{\nu}+Pg^{\mu\nu},$ (3)

where the specific enthalpy is defined by

$h=1+\epsilon+\frac{P}{\rho_{0}}.$ (4)

The active mass density is then given by

$\rho_{\rm{active}}=T-2T^{0}_{0}=T^{i}_{i}-T^{0}_{0}=\rho_{0}h\frac{1+v^{2}}{1-v^{2}}+2P.$
(5)

In Eqs. (3), (4) and (5), $\rho_{0}$ is the rest mass density which is
proportional to the number density of baryons of the fluid, $P$ is the
pressure, $u^{\mu}$ is the four velocity of a fluid element with respect to an
Eulerian observer, $\epsilon$ is the specific internal energy, $T$ is the
trace of the energy-momentum tensor ($T=g_{\mu\nu}T^{\mu\nu}$), and $v$ is the
three dimensional fluid velocity. Unlike the Newtonian case,
$\rho_{\rm{active}}$ includes all sources of energy.

The equations governing the motion of the fluid matter can be derived from the
conservations laws for the energy-momentum tensor and the fluid’s matter
current, i.e., $\nabla_{\mu}T^{\mu\nu}=0$ and $\nabla_{\mu}J^{\mu}=0$. In the
ADM decomposition of spacetime (Arnowitt et al., 1962), the metric
($g_{\mu\nu}$) can be expressed in the following form by considering the
foliation of spacetime using 3-dimensional hypersurfaces defined by $t={\rm
const.}$:

$ds^{2}=-\alpha^{2}dt^{2}+\gamma_{ij}\left(dx^{i}+\beta^{i}dt\right)\left(dx^{j}+\beta^{j}dt\right),$
(6)

Here $\gamma_{ij}$ is the spatial metric, defined on each hypersurface, while
$\alpha$ and $\beta^{i}$ are known as the lapse, and shift vector,
respectively, and encode the 4-fold coordinate freedom of general relativity.

Flux-conservative formulations of hydrodynamics have been applied very
successfully in computational fluid dynamics. To cast the fluid equations in
flux-conservative form we first define so-called conservative variables ($q$)
in terms of the original hydrodynamic variables (so-called primitive
variables, $w$),

$q=\left(\begin{array}[]{c}D\\\ S_{i}\\\
\tau\end{array}\right)=\left(\begin{array}[]{c}\rho_{0}W\\\
\rho_{0}hW^{2}v_{i}\\\
\rho_{0}hW^{2}-P-D\end{array}\right),\,w=\left(\begin{array}[]{c}\rho_{0}\\\
v^{i}\\\ P\end{array}\right),$ (7)

where $W=1/\sqrt{1-\gamma_{ij}v^{i}v^{j}}$. With these definitions, and with
the metric (Eq. (6)), we can then write the Euler equation as (Font et al.,
2000)

$\frac{\partial\left(\sqrt{\gamma}q\right)}{\partial
t}+\frac{\partial\left(\sqrt{-g}f^{i}\right)}{\partial
x^{i}}=\sqrt{-g}\Sigma,$ (8)

where the fluxes $f^{i}$ and the sources $\Sigma$ are given by

$\displaystyle f^{i}$ $\displaystyle=$
$\displaystyle\left[\begin{array}[]{ccccc}D\left(v^{i}-\frac{\beta^{i}}{\alpha}\right)\\\
S_{j}\left(v^{i}-\frac{\beta^{i}}{\alpha}\right)+P\delta^{i}_{j}\\\
\tau\left(v^{i}-\frac{\beta^{i}}{\alpha}\right)+Pv^{i}\end{array}\right],$
(12) $\displaystyle\Sigma$ $\displaystyle=$
$\displaystyle\left[\begin{array}[]{ccccc}0\\\
T^{\mu\nu}\left(\partial_{\mu}g_{\mu j}-\Gamma^{\lambda}_{\mu\nu}g_{\lambda
j}\right)\\\ \alpha\left(T^{\mu
0}\partial_{\mu}\left(\ln\alpha\right)-\Gamma^{0}_{\mu\nu}T^{\mu\nu}\right)\end{array}\right].$
(16)

Here $\sqrt{\gamma}$ and $\sqrt{-g}$ are the determinants of $\gamma_{ij}$ and
$g_{\mu\nu}$, respectively, and are related by
$\sqrt{-g}=\alpha\sqrt{\gamma}$. It is well known that for a perfect fluid,
the system of equations derived from the conservation laws is not closed: the
number of dynamical equations is always less than the number of unknowns.

As is also well known, the equation of state (hereafter EOS) for the fluid
provides an additional equation, but in the general case it also introduces
other unknowns. In order to completely close the hydrodynamical equations, an
energy balance equation is often used. However, under certain circumstances,
we can adopt rather simple EOSs that do not introduce any further variables:
adiabatic and isothermal EOSs provide specific examples.

Realistic EOSs are usually determined by theoretical calculations and
experimental measurements. However, there are physical regimes where our
understanding of the nature of the matter is quite incomplete. Specifically,
in the case where the matter density is significantly above nucleon density,
there remain large uncertainties in the correct EOS. Thus, for example, the
EOS at the core of neutron stars is still not very well understood. Here we
ignore these difficulties, and for the purpose of testing our code, use two
types of very simple EOS. The first is the ideal gas EOS which can be written
in the following form:

$P=\left(\Gamma-1\right)\rho_{0}\epsilon,$ (17)

and corresponds to the isothermal EOS. We use this EOS in the shock tube tests
described in (section 4). The second EOS results from the isentropic
assumption, whereby Eq. (17) becomes the polytropic EOS:

$P=K\rho_{0}^{1+\frac{1}{N}}.$ (18)

Here $K$ and $N$ are the polytropic constant and index respectively. The
polytropic EOS of state is the generalized form of the adiabatic one; a fluid
which is governed by it does not generate entropy, and shock formation is thus
generically prohibited. We use this EOS in the pulsation mode test (section 6
and Appendix A).

Using the above formulation, we are now ready to describe in detail the
pseudo-Newtonian hydrodynamical equations used in our code. We limit our study
here to axisymmetric systems, and adopt cylindrical coordinates ($R,Z,\phi$)
such that

$ds^{2}=-(1+2\Phi)dt^{2}+\frac{1}{1+2\Phi}\left(dR^{2}+dZ^{2}+R^{2}d\phi^{2}\right).$
(19)

The lapse function and shift vector are thus given by $\alpha=\sqrt{1+2\Phi}$
and $\beta^{i}=0$. In addition, we enforce the equatorial symmetry at $z=0$
since the phenomena involving $l=\rm{odd}$ modes are not dominant in most
cases of neutron star dynamics, where $l$ is from the spherical harmonics. In
this coordinate system, the conservative and primitive variables are

$q=\left(\begin{array}[]{c}D\\\ S_{R}\\\ S_{Z}\\\ S_{\phi}\\\
\tau\end{array}\right)=\left(\begin{array}[]{c}\rho_{0}W\\\
\rho_{0}hW^{2}v_{R}\\\ \rho_{0}hW^{2}v_{Z}\\\ \rho_{0}hW^{2}v_{\phi}\\\
\rho_{0}hW^{2}-P-D\end{array}\right),\,w=\left(\begin{array}[]{c}\rho_{0}\\\
v^{R}\\\ v^{Z}\\\ v^{\phi}\\\ P\end{array}\right).$ (20)

The final form of the hydrodynamical equations then becomes

$\frac{\partial\left(\sqrt{\gamma}q\right)}{\partial
t}+\frac{\partial\left(\sqrt{-g}f^{R}\right)}{\partial
R}+\frac{\partial\left(\sqrt{-g}f^{Z}\right)}{\partial Z}=\sqrt{-g}\Sigma$
(21)

where

$\displaystyle f^{R}$ $\displaystyle=$
$\displaystyle\left[\begin{array}[]{ccccc}Dv^{R}\\\ S_{R}v^{R}+P\\\
S_{Z}v^{R}\\\ S_{\phi}v^{R}\\\ \tau v^{R}+Pv^{R}\end{array}\right],$ (27)
$\displaystyle f^{Z}$ $\displaystyle=$
$\displaystyle\left[\begin{array}[]{ccccc}Dv^{Z}\\\ S_{R}v^{Z}\\\
S_{Z}v^{Z}+P\\\ S_{\phi}v^{Z}\\\ \tau v^{Z}+Pv^{Z}\end{array}\right],$ (33)
$\displaystyle\Sigma$ $\displaystyle=$
$\displaystyle\left[\begin{array}[]{ccccc}0\\\
-\frac{\rho_{\rm{active}}}{1+2\Phi}\frac{\partial\Phi}{\partial
R}+\frac{S_{\phi}v^{\phi}}{R}+\frac{P}{R}\\\
-\frac{\rho_{\rm{active}}}{1+2\Phi}\frac{\partial\Phi}{\partial Z}\\\ 0\\\
-\left(S_{R}\frac{\partial\Phi}{\partial R}+S_{Z}\frac{\partial\Phi}{\partial
Z}\right)\end{array}\right].$ (39)

Using Eq. (19) we have $\sqrt{\gamma}=R\left(1+2\Phi\right)^{-3/2}$, and
$\sqrt{g}=R\left(1+2\Phi\right)^{-1}$. In obtaining the expressions in Eq.
(39) we have used the assumption of slow changes of the potential relative to
the gradients ($\frac{\partial\Phi}{\partial t}\ll\frac{\partial\Phi}{\partial
R}$ or $\frac{\partial\Phi}{\partial Z}$). Recently, Nagakura et al. (2011)
used a similar method in the context of jet propagation in a uniform medium,
but adopted a slightly different linear momentum equation than ours. (See Eq.
(39) and compare with Eqs. (2) and (3) in Nagakura et al., 2011).

Finally, the gravitational Poisson equation in our coordinate system is

$\frac{1}{R}\frac{\partial}{\partial R}\left(R\frac{\partial\Phi}{\partial
R}\right)+\frac{\partial^{2}\Phi}{\partial Z^{2}}=4\pi\rho_{\rm{active}}.$
(40)

Note that the second component of $\Sigma$ contains terms which, individually,
become singular on the axis of symmetry ($R=0$). In addition, there are other
terms in the equations of motion that need to be treated carefully as $R\to
0$. This is done by demanding regularity at the axis, and by considering the
parity of each function, with respect to $R$, in that limit. In particular,
$\rho_{0}$, $v^{Z}$, $v^{\phi}$, $P$, $D$, $S_{Z}$, $S_{\phi}$, and $\tau$ are
all even functions of $R$ as $R\to 0$, while $v^{R}$ and $S_{R}$ are odd.
Taking this into account, Eq. (21) and $\Sigma$ in Eq. (39) become

$\frac{\partial\left(\sqrt{\gamma^{\prime}}q\right)}{\partial
t}+2\frac{\partial\left(\sqrt{-g^{\prime}}f^{R}\right)}{\partial
R}+\frac{\partial\left(\sqrt{-g^{\prime}}f^{Z}\right)}{\partial
Z}=\sqrt{-g^{\prime}}\Sigma$ (41)

and

$\Sigma=\left[\begin{array}[]{ccccc}0\\\ 0\\\
-\frac{\rho_{\rm{active}}}{1+2\Phi}\frac{\partial\Phi}{\partial Z}\\\ 0\\\
-S_{Z}\frac{\partial\Phi}{\partial Z}\end{array}\right],$ (42)

where $\sqrt{\gamma^{\prime}}=\left(1+2\Phi\right)^{-3/2}$ and
$\sqrt{-g^{\prime}}=\left(1+2\Phi\right)^{-1}$. The coefficient of
$\frac{\partial\left(\sqrt{-g}f^{R}\right)}{\partial R}$ in Eq. (41) becomes
$2$ instead of $1$, while the other variables, such as $q$, $f^{R}$ and
$f^{Z}$, are unchanged from Eq. (39).

Finally, using L’Hopital’s theorem at $R=0$, the singular term
$R^{-1}\partial\Phi/\partial R$ in the Poisson equation (40) is replaced by
$\partial^{2}\Phi/\partial R^{2}$.

##  3 Numerical Methods

In this section, we describe our numerical methods for solving the coupled
hydrodynamical and Poisson equations. We mainly use the finite volume methods
for the hydrodynamical equations and the a finite difference approach for the
Poisson equation. In the finite volume method, each grid cell represents
volume averaged hydrodynamic quantities i.e., $\bar{q}=\frac{1}{\Delta V}\int
qdV$. After applying the finite volume method our hydro equations can be
reduced to Riemann problems which consider the time evolution of initial
conditions given by two distinct states that adjoin at some interface (so that
there are, in general, discontinuities across one or more physical quantities
at the interface). A very important property of the finite volume method is
that it maintains the local conservation properties of the flow in the
computational grid.

In the dynamics of compressible fluids, we inevitably encounter discontinuous
behaviors such as shocks, rarefactions or contact discontinuities. To treat
such discontinuities without introducing numerical instabilities or spurious
oscillations, we use High Resolution Shock Capturing (HRSC) techniques that
generically reduce the order of accuracy of the numerical scheme near
discontinuities or when one or more of the fluid variables are at a local
maximum. A key ingredient to the success of the HRSC methods is the
calculation of fluxes through cell boundaries. To compute these fluxes we need
approximate values for the primitive variables at the cell boundaries. We have
implemented second order slope limiters such as minmod (van Leer, 1979),
monotonized central difference (MC hereafter, van Leer, 1977) and superbee
(Roe, 1985), as well as a third order slope limiter proposed by Shibata (2003)
and which is based on the minmod function (3minmod hereafter). Other
reconstruction methods such as the third order Piecewise Parabolic Method (PPM
hereafter, Colella & Woodward, 1984), Essentially Non-Oscillatory method (ENO,
Harten et al., 1987) and Weighted ENO (WENO, Liu et al., 1994; Jiang & Shu,
1996), which has an arbitrary order of accuracy, were also implemented.

In the implementation of HRSC schemes it is not efficient to exactly solve the
Riemann problems which arise since an excessively large amount of
computational resources per cell are then needed to calculate the fluxes.
Thus, an approximate calculation of fluxes is performed. We implemented the
following three schemes: Roe (Roe, 1981), Marquina (Donat & Marquina, 1996;
Donat, 1998), and HLLE (Harten et al., 1983; Einfeldt, 1988; Einfeldt et al.,
1991) approximations. The Roe approximation is based on the Rankine-Hugoniot
jump condition and Marquina’s approach generalizes Roe’s scheme. The HLLE
algorithm comes from a very simple two wave approximation and produces the
most dissipative and stable results. We have mainly used HLLE for reducing
computational cost but found that our results were not significantly
influenced by the type of flux approximation used.

In order to solve the Poisson equation (which is elliptic) for the
gravitational potential we use the multigrid method, which can quickly reduce
low frequency error components in the solution by adopting hierarchical grid
levels (Brandt, 1977). One of the difficulties we often encounter with the
Poisson equation is in the proper implementation of the boundary conditions.
For example, one of the natural boundary conditions is $\Phi=0$ at $\infty$,
but in the coordinates adopted in the previous section the computational
domain cannot reach spatial infinity. We thus now refer to our previous
coordinates as $(r,z)$ and introduce new coordinates $(R,Z)$ which compactify
the spatial domain, mapping the infinities in each spatial direction to finite
coordinate values. Specifically, we choose the same type of compactification
for both $r$ and $z$ coordinates, namely a tangent function, but allow a
certain portion of the domain to remain ”’uncompactified”:

$\displaystyle r=\left\\{\begin{array}[]{ll}R&\mbox{if $R\leq r_{0}$}\\\
r_{0}+r_{1}\tan\left(\frac{R-r_{0}}{r_{1}}\right)&\mbox{if
$R>r_{0}$}\end{array}\right.,$ (45) $\displaystyle
z=\left\\{\begin{array}[]{ll}Z&\mbox{if $Z\leq z_{0}$}\\\
z_{0}+z_{1}\tan\left(\frac{Z-z_{0}}{z_{1}}\right)&\mbox{if
$Z>z_{0}$}\end{array}\right..$ (48)

Here, the four parameters $z_{0}$, $z_{1}$, $r_{0}$ and $r_{1}$ control the
compactification, and we chose this specific form for the coordinate
transformation since it guarantees that the compactified coordinates smoothly
transition to the original ones near the origin. We note that we solve the
hydrodynamical and gravitational equations on separate spatial domains:
$[0:r_{0},0:z_{0}]$ for the hydrodynamic calculations and
$[0:r_{0}+\frac{2}{\pi}r_{1},0:z_{0}+\frac{2}{\pi}z_{1}]$ for the computation
of the gravitational potential. There ranges correspond to $[0:r_{0},0:z_{0}]$
and $[0:\infty,0:\infty]$, respectively, in the original cylindrical
coordinates $(r,z)$. In the compactified coordinates, the Poisson equation is
written as

$\displaystyle\frac{1}{rf(R)}\frac{\partial}{\partial
R}\left(\frac{r}{f(R)}\frac{\partial\Phi(R,Z)}{\partial R}\right)$
$\displaystyle+\frac{1}{g(Z)}\frac{\partial}{\partial
Z}\left(\frac{1}{g(Z)}\frac{\partial\Phi(R,Z)}{\partial
Z}\right)=4\pi\rho_{\rm{active}},$ (49)

where $f(R)$ and $g(Z)$ are given by

$\displaystyle f(R)=\frac{dr}{dR}=\left\\{\begin{array}[]{ll}1&\mbox{if $R\leq
r_{0}$}\\\ \sec^{2}\left(\frac{R-r_{0}}{r_{1}}\right)&\mbox{if
$R>r_{0}$},\end{array}\right.$ (52) $\displaystyle
g(Z)=\frac{dz}{dZ}=\left\\{\begin{array}[]{ll}1&\mbox{if $Z\leq z_{0}$}\\\
\sec^{2}\left(\frac{Z-z_{0}}{z_{1}}\right)&\mbox{if
$Z>z_{0}$}\end{array}\right..$ (55)

As just noted, the domain for the hydrodynamical calculation is finite, i.e.
we do not solve the hydrodynamical equations on the full compactified domain,
and we thus must be careful to choose values of $r_{0}$ and $z_{0}$ large
enough so that there is no outflux of matter through the $r=r_{0}$ and/or
$z=z_{0}$ boundaries. In our code we set $r_{0}=z_{0}=\eta r_{e}$ where $\eta$
is a free parameter and $r_{e}$ is the equatorial radius of the rotating star
as obtained from the procedure we use to calculate the initial stellar model.
For the pulsation mode test described in section 6, a typical choice is
$\eta=2$. This means that the hydrodynamical computational domain extends
twice the distance of the stellar radius in both the $R$ and $Z$ directions:
this choice is found to be sufficient for our study. The values of $r_{1}$ and
$z_{1}$ are automatically determined by requiring the multigrid domain to be 2
times larger than the size of hydrodynamic domain in compactified coordinates,
i.e., $r_{1}=\frac{\pi}{2}r_{0}$ and $z_{1}=\frac{\pi}{2}z_{0}$.

In our multigrid algorithm, we use line relaxation for our basic smoother,
whereby all grid point values given by $R={\rm const.}$ or $Z={\rm const.}$
are updated simultaneously (constant-$R$ and constant-$Z$ sweeps are
alternated). We cannot use point-wise relaxation since, as is well known, such
a technique is not a good smoother when there is significant anisotropy in the
coefficients of the second derivative terms in the elliptic operator being
treated. This is the case in our compactified coordinate system, particularly
near the domain boundaries. We have used second- and fourth-order finite-
difference approximations to the Poisson equation, and these lead to tri- and
pentadiagonal linear systems, respectively, that must be solved to implement
the line relaxations. We use the routines DGTSV (tridiagonal) and DGBSV
(banded/pentadiagonal) routines from LAPACK to perform these solutions.

In order to integrate the discretized hydrodynamical equations, we use the
method of lines (MOL), transforming our partial differential equations in time
and space to ordinary differential equations (ODEs) with respect to the time.
To solve these ODEs, we then employ second and third order Runge-Kutta
methods, which are known to have the Total Variation Diminishing (TVD)
property.

##  4 Shock tube Tests

In order to verify the accuracy and the convergence of our numerical code, we
first carry out rigorous tests using initial configurations having analytic
solutions. In this section, we present the results of such tests for the case
where there is no self-gravity (i.e. pure hydrodynamics). Another test of the
entire code—including our treatment of the gravitational field—is described in
the next section.

Shock tube tests are Riemann problems where the initial configuration of the
fluid is given by two states having, in general, different densities,
pressures and velocities, on the left and right halves of the tube. Three
possible distinct features emerge from the subsequent evolution: a shock, a
rarefaction fan, and a contact discontinuity. We carried out 1D and 2D
numerical simulations with 3 different parameter sets previously used by Zhang
& MacFadyen (2006). These parameters are listed in Table 1, where superscripts
$R$ and $L$ represent the fluid states in the right and left halves,
respectively, of the tube.

Table 1: Initial values of physical quantities for Shock tube tests (Riemann problem). Problem | $\Gamma$ | $\rho_{0}^{L}$ | $\rho_{0}^{R}$ | $v^{L}$ | $v^{R}$ | $P^{L}$ | $P^{R}$
---|---|---|---|---|---|---|---
1 | 5/3 | $10.0$ | $1.0$ | $0.0$ | $0.0$ | $13.33$ | $10^{-8}$
2 | 5/3 | $1.0$ | $1.0$ | $0.0$ | $0.0$ | $1000.0$ | $10^{-2}$
3 | 4/3 | $1.0$ | $1.0$ | $0.9$ | $0.0$ | $1.0$ | $10.0$

###  4.1 1D test in Cartesian Coordinates

We first we present the results of our 1D tests. In these tests, we carefully
examine how the distinct features predicted by the analytic solutions are
reproduced by different methods, and measure the accuracy and convergence rate
of the various solutions obtained.

In problem 1, the initial discontinuity gives 3 different types of solutions
(shock, rarefaction, contact discontinuity). Figure 1 shows the results at
$t=0.4$ obtained using four different methods of reconstruction: minmod (top
left), MC (top right), 3minmod (bottom left) and PPM (bottom right). We
observe that the minmod method is quite dissipative, yielding rather smooth
solutions that cannot accurately describe the shockwave. We also find that at
low resolution the height of the shock is not well reproduced if we use the
minmod method. MC and 3minmod give almost similar results, while PPM shows the
best behaviour near the shock.

00.20.40.60.8100.20.40.60.81 $\rho_{0}$, $P$, $v$ $x$-axis$\rho_{0}/10$$P/15$$v$ |  00.20.40.60.8100.20.40.60.81 $\rho_{0}$, $P$, $v$ $x$-axis$\rho_{0}/10$$P/15$$v$
---|---
00.20.40.60.8100.20.40.60.81 $\rho_{0}$, $P$, $v$ $x$-axis$\rho_{0}/10$$P/15$$v$ |  00.20.40.60.8100.20.40.60.81 $\rho_{0}$, $P$, $v$ $x$-axis$\rho_{0}/10$$P/15$$v$
Figure 1: One dimensional shock tube test of problem 1 at $t=0.4$ with
different reconstruction methods: minmod(top-left), MC (top right), 3minmod
(bottom left), and PPM (bottom right). The initial discontinuity is at
$x=0.5$. We use 512 uniform grid points. The numerical results are shown in 3
different colours: rest mass density (pink), pressure (red) and velocity
(blue). The solid lines show the analytic solutions.

The second test problem (Problem 2) is the so-called blast wave test which
produces a very sharp and thin shell in density between the shock and contact
discontinuity. Generally, numerical codes are not able to perfectly resolve
this very thin shell because it can span only a few grid cells, even in very
high resolution calculations. Nonetheless, this test provides insight as to
how well a code can handle such a feature. As can see in Figure 2, PPM again
gives the best results. although it still shows large errors at the shock.

00.20.40.60.8100.20.40.60.81 $\rho_{0}$, $P$, $v$ $x$-axis$\rho_{0}/15$$P/1000$$v$ |  00.20.40.60.8100.20.40.60.81 $\rho_{0}$, $P$, $v$ $x$-axis$\rho_{0}/15$$P/1000$$v$
---|---
00.20.40.60.8100.20.40.60.81 $\rho_{0}$, $P$, $v$ $x$-axis$\rho_{0}/15$$P/1000$$v$ |  00.20.40.60.8100.20.40.60.81 $\rho_{0}$, $P$, $v$ $x$-axis$\rho_{0}/15$$P/1000$$v$
Figure 2: Same as Fig. 1 for problem 2.

The third problem generates a strong reverse shock but numerical solution has
oscillatory features near the shock front. Generally speaking, the oscillation
can be easily damped out if the numerical scheme is significantly dissipative.
Numerical dissipation also tends to weaken the sharpness of the discontinuity.
In Figure 3, one can see that the minmod methods, which, as already noted, is
the most dissipative of the techniques we use, gives relatively small
amplitude oscillations, except near the discontinuity. The more non-
dissipative methods describes the shock features well, but produce rather
large amplitude oscillatory behavior.

00.20.40.60.8100.20.40.60.81 $\rho_{0}$, $P$, $v$ $x$-axis$\rho_{0}/10$$P/25$$v$ |  00.20.40.60.8100.20.40.60.81 $\rho_{0}$, $P$, $v$ $x$-axis$\rho_{0}/10$$P/25$$v$
---|---
00.20.40.60.8100.20.40.60.81 $\rho_{0}$, $P$, $v$ $x$-axis$\rho_{0}/10$$P/25$$v$ |  00.20.40.60.8100.20.40.60.81 $\rho_{0}$, $P$, $v$ $x$-axis$\rho_{0}/10$$P/25$$v$
Figure 3: Same as Figs. 1 and 2 for problem 3.

To quantify the deviation of our numerical results from the analytic
solutions, we use the $L_{1}$ norm of the errors, defined by
$L_{1}=\sum_{i=1}^{N}\Delta x_{i}|q_{i}-q(x_{i})|$, where $q(x_{i})$ is the
value of the analytic solution at point $x_{i}$. We summarize the results in
Table 2. The convergence rate ($\log_{2}\left[L_{1}^{2h}/L_{1}^{h}\right]$) in
the table should be close to 1, which corresponds to the $1^{\textrm{st}}$
order nature of the HRSC scheme near the shock where the most of the
$L_{1}$-norm error occurs. However, it can deviate from that value due to the
oscillatory features near the shock.

Table 2: The $L_{1}$ norm of the error and its convergence rate for each of the test problems using different resolutions and different reconstruction schemes. | N
---|---
| 64 | 128 | 256 | 512 | 1024 | 2048
Problem 1 | minmod |  $L_{1}$ norm ($\times 10^{-2}$) | $25.7$ | $16.5$ | $9.48$ | $5.02$ | $2.66$ | $1.49$
convergence rate | - | $0.64$ | $0.80$ | $0.92$ | $0.91$ | $0.84$
MC | - | $15.3$ | $9.43$ | $5.25$ | $2.79$ | $1.46$ | $0.830$
- | - | $0.70$ | $0.85$ | $0.91$ | $0.93$ | $0.82$
3minmod | - | $17.8$ | $11.0$ | $5.82$ | $2.99$ | $1.51$ | $0.816$
- | - | $0.69$ | $0.92$ | $0.96$ | $0.98$ | $0.89$
PPM | - | $12.3$ | $6.55$ | $3.43$ | $1.74$ | $0.877$ | $0.431$
- | - | $0.91$ | $0.93$ | $0.98$ | $0.99$ | $1.0$
Problem 2 | - | - | $30.1$ | $21.0$ | $20.1$ | $15.8$ | $10.9$ | $6.93$
- | - | $0.52$ | $0.061$ | $0.34$ | $0.54$ | $0.65$
- | - | $27.8$ | $18.2$ | $14.9$ | $10.4$ | $6.28$ | $3.77$
- | - | $0.61$ | $0.29$ | $0.52$ | $0.73$ | $0.74$
- | - | $28.3$ | $17.9$ | $13.7$ | $8.84$ | $5.05$ | $2.72$
- | - | $0.66$ | $0.39$ | $0.63$ | $0.81$ | $0.89$
- | - | $29.5$ | $17.9$ | $12.7$ | $7.79$ | $3.73$ | $2.13$
- | - | $0.73$ | $0.49$ | $0.71$ | $1.1$ | $0.81$
Problem 3 | - | - | $15.5$ | $10.0$ | $6.19$ | $3.65$ | $2.37$ | $1.58$
- | - | $0.63$ | $0.69$ | $0.76$ | $0.63$ | $0.59$
- | - | $14.9$ | $7.73$ | $5.40$ | $2.72$ | $1.64$ | $1.04$
- | - | $0.95$ | $0.52$ | $0.99$ | $0.73$ | $0.66$
- | - | $13.0$ | $6.97$ | $4.40$ | $2.25$ | $1.35$ | $0.867$
- | - | $0.90$ | $0.66$ | $0.97$ | $0.74$ | $0.63$
- | - | $7.26$ | $3.93$ | $2.41$ | $1.08$ | $0.547$ | $0.393$
- | - | $0.89$ | $0.70$ | $1.16$ | $0.98$ | $0.47$
0.010.11PPM3-minmodMCminmod $\rm{L}_{1}\rm{errornorm}$ Problem 1Problem
2Problem 3
---
0.40.60.811.21.4PPM3-minmodMCminmod Convergence Rate Problem 1Problem 2Problem
3
Figure 4: The $L_{1}$ norm (top) and convergence rate (bottom) when the number
of grid points is 512 with different reconstruction methods. Three different
shock tube problems are shown with different colors (problem 1: blue, problem
2: red and problem 3: sky blue).

Figure 4 shows the $L_{1}$ norms and convergence rates for each problem when
the grid resolution is $\rm{N}=512$. Although no single method stands out in
our 1D shock tube tests, we conclude from from the values of the $L_{1}$ error
norms and convergence rates that PPM gives the most promising results.

###  4.2 2D test in Cylindrical Coordinates

Since the cylindrical coordinate system we have adopted is curvilinear,
1-dimensional shock tube tests are not sufficient for assessing our code’s
accuracy and convergence. In Cartesian coordinates, fluxes between cells which
have the same state cancel out. For example, if we carry out the shock tube
test in the $x$-direction, then the fluxes in the $y$ and $z$ directions are
identical in every grid cell, meaning that the net flux is 0. Therefore, 1D
shock tube tests performed with codes that use 2- or 3D Cartesian coordinates
produce exactly the same results as a 1D code. However, in cylindrical
coordinates, fluxes do not cancel in this way, but rather are balanced by
source terms. This difference may give additional non-physical effects,
especially near discontinuities.

Therefore, we carried out the first of the shock tube tests listed in Table 1
in cylindrical coordinates, where we placed the discontinuity on the Z=0
plane. Figure 5 shows the solution resulting solution on the $Z$-axis. If we
use the minmod method, the 2D results are similar to the 1 dimensional ones.
In addition, although PPM produces better results minmod, it cannot produce
the sharp features of the shock seen in the 1D test: this is due to the
dissipation caused by the imbalance between the net flux and the source term.
We can also see that 3minmod and PPM yield quite similar results. We checked
the differences in solutions at different $z={\rm const.}$ planes and found
that they are negligibly small ($\sim 10^{-13}$) compared to the truncation
errors. Overall, however, although the 2D results show more dissipation than
the 1D ones, the relative differences in the solutions are not significant
(for our purposes). In particular, both agree acceptably with the analytic
forms.

00.20.40.60.8100.20.40.60.81

$\rho_{0}/10$

$x$-axisminmod3minmodPPM

Figure 5: The solution on the axis at t=0.4 in problem 1 in section 4.1. The 3
different reconstruction methods (minmod: blue, 3-minmod: sky blue and PPM:
magenta) are shown.

##  5 Stationary Star Test

The tests just reported did not involve the effects of the gravitational
field. In this section, we test our treatment of the Poisson equation for the
gravitational potential as well as the hydrodynamics.

With an ideal code, the evolution of a stationary star should also be
stationary. However, in practice, all codes that dynamically evolve stationary
states show some level of fluctuation due to finite grid resolution and
intrinsic errors in the numerical scheme used. In this section, we show the
time evolution of the physical quantities of non-rotating and the rotating
stars, and investigate the dependence of this time behaviour by changing the
resolution of the simulations. Specifically, we use 3 different grid
resolutions : $65\times 65$, $129\times 129$ and $257\times 257$, where $1/2$
of the grid points span the star at the equator.

Our initial models of rotating stars are generated using Hachisu’s Self-
Consistent Field (HSCF: Hachisu, 1986a, b) method—details of the procedure are
described in Kim et al. (2009). In order to generate equilibrium models, we
choose 1) the maximum rest mass density, $\rho_{0}^{\rm{max}}$) 2) the
rotation parameter, $A$, which describes the differential rotation and 3) the
axis ratio which determines how fast the star is rotating. We must also
specify the equation of state (EOS) in our construction of the initial model.
Here, we used the polytropic EOS Eq. (18) with $K=100$ and $N=1$. We choose a
maximum density value of $\rho_{0}^{\rm{max}}=1.28\times 10^{-3}$, which, with
this EOS, produces a 1.4$\rm{M}_{\odot}$star in the non-rotating case. For the
rotating models, we only consider rigid body rotation, which is obtained when
we chose a very large value of $A$. The axis ratio is specified to be $0.75$
resulting in an orbital frequency of $611\rm{Hz}$.

Even with our use of the multigrid technique—which is generally an efficient
method for solving elliptic equations—we still find solution of the Poisson
equation for the gravitational potential to be computationally expensive. We
thus calculate $\Phi$ only every $40$ time steps to reduce the time spent in
the Poisson solver, and find that this produces results which are nearly
equivalent to those obtained when the Poisson equation is solved at each time
step. However, we use time-extrapolated values for the gravitational potential
at the time steps between solves of the Poisson equation in order to avoid
discontinuities in the primitive variables, when abrupt changes of the
gravitational potential occur. We find that these discontinuities give rise to
very unnatural dissipative effects in the simulation resulting, for example,
in a rapid decay in the amplitude of radial oscillations, even when radial
perturbations are explicitly introduced.

-0.01-0.00500.0050.010246810 $\left[\rho_{0}^{\rm{max}}(t)-\rho_{0}^{\rm{max}}(0)\right]/\rho_{0}^{\rm{max}}(0)$ time(ms)$\mathcal{R}=3.05\times 10^{-7}$$65\times 65$$129\times 129$$257\times 257$ |  -0.01-0.00500.0050.010246810 $\left[\rho_{0}^{\rm{max}}(t)-\rho_{0}^{\rm{max}}(0)\right]/\rho_{0}^{\rm{max}}(0)$ time(ms)$\mathcal{R}=3.95\times 10^{-6}$$65\times 65$$129\times 129$$257\times 257$
---|---
-0.01-0.00500.0050.010246810 $\left[\rho_{0}^{\rm{max}}(t)-\rho_{0}^{\rm{max}}(0)\right]/\rho_{0}^{\rm{max}}(0)$ time(ms)$\mathcal{R}=1.94\times 10^{-7}$$65\times 65$$129\times 129$$257\times 257$ |  -0.01-0.00500.0050.010246810 $\left[\rho_{0}^{\rm{max}}(t)-\rho_{0}^{\rm{max}}(0)\right]/\rho_{0}^{\rm{max}}(0)$ time(ms)$\mathcal{R}=3.42\times 10^{-6}$$65\times 65$$129\times 129$$257\times 257$
Figure 6: The time evolution of the maximum rest mass density
changes($\left[\rho_{0}^{\rm{max}}(t)-\rho_{0}^{\rm{max}}(t=0)\right]/\rho_{0}^{\rm{max}}(t=0)$)
with different resolutions ($65\times 65$: red, $129\times 129$: dark blue and
$257\times 257$: sky blue). Top figure shows results when we fix the metric
(Cowling approximation) while the bottom one shows the case where we consider
the fully coupled dynamics. In the left panel, we show the figures for a
spherical (non-rotating) star while the right panel shows the corresponding
figures for a rigidly rotating star with axis ratio$=0.75$, which give a
rotational frequency of $611\rm{Hz}$

Figure 6 shows the time evolution of the relative changes of the maximum
density
($\left[\rho_{0}^{\rm{max}}(t)-\rho_{0}^{\rm{max}}(0)\right]/\rho_{0}^{\rm{max}}(0)$)
for non-rotating (left panel) and rigidly rotating (right panel) stars. For
the stationary stars, we use the Cowling approximation, which assumes the
gravitational potential is fixed. This gives efficient evolution of the stars,
and can also be used as a testbed for fully coupled evolutions. The results
computed using the Cowling approximation are shown in the top figures. The
maximum density slowly increases with time for the rotating star while it
decreases for the non-rotating star. For grid resolutions greater than
$65\times 65$ the rate of change is almost independent of resolution for the
spherical star, but a slow decrease with resolution is seen for the rotating
star. We define the following dimensionless rate of change:

$\mathcal{R}=\left|t_{\rm{dyn}}\frac{d\ln\rho_{0}^{\rm{max}}}{dt}\right|,$
(56)

where we use $t_{\rm{dyn}}=1/\sqrt{\rho_{0}^{\rm{max}}}$ for simplicity. We
use this quantity—as computed from the highest resolution simulations—as a
label in the figures. The values of $\mathcal{R}$ are within $3\times 10^{-7}$
for non-rotating star and are about 10 times larger for the rotating star,
again with a maximum resolution of $257\times 257$. The inverse of
$\mathcal{R}$ can be interpreted as the time (in units of the dynamical time)
that the simulation could be carried out until the results deviate from the
true solution by $O(1)$. Our results indicate that the error would become
$\sim$ 1% in 30,000 and 3,000 dynamical times for non-rotating and rotating
stars, respectively. We also carried out very long time simulations and found
that $\mathcal{R}$ becomes smaller even though it appears to be almost
constant in the figures. From these results, we conclude that we can use the
code to evolve stellar configurations for several thousand or more dynamical
times.

It is also very important to check the constancy of the conserved quantities
with respect to simulation time. In our formulation we have 2 conserved
quantities: the total rest mass, $M_{0}$, and the total angular momentum, $J$,
which are computing using

$\displaystyle M_{0}=\int
DdV^{(3)}=2\pi\int\frac{\rho_{0}W}{\left(1+2\Phi\right)^{3/2}}RdRdZ,$ (57)
$\displaystyle J=\int
S_{\phi}dV^{(3)}=2\pi\int\frac{\rho_{0}hW^{2}v^{\phi}}{\left(1+2\Phi\right)^{5/2}}R^{3}dRdZ,$
(58)

respectively, and where $dV^{(3)}$ denotes the 3-dimensional volume element.

-0.0001-5e-0505e-050.00010246810 $\left[M_{0}(t)-M_{0}(0)\right]/M_{0}(0)$ time(ms)$\mathcal{R}_{M}=4.79\times 10^{-7}$$65\times 65$$129\times 129$$257\times 257$
---
-0.001-0.0008-0.0006-0.0004-0.000200246810 $\left[J(t)-J(0)\right]/J(0)$ time(ms)$\mathcal{R}_{J}=1.67\times 10^{-6}$$65\times 65$$129\times 129$$257\times 257$
Figure 7: The deviation of total rest mass (upper panel) and angular momentum
(lower panel), which should remain constant, from their initial values with
time, for the same models shown in Fig. 6. Results computed with three
different resolutions ($65\times 65$:red, $129\times 129$:dark blue and
$257\times 257$:sky blue) are shown.

Figure 7 shows the time evolution of these two conserved quantities: total
rest mass (upper panel) and total angular momentum (lower panel). We show the
results only from the rotating star since there is, of course, no angular
momentum for non-rotating stars. The deviation of the total rest mass from the
initial value has two features: short-term fluctuations and long term average
behavior. The shot-term fluctuations depend on the grid resolution, but the
average slopes are almost independent of the resolution. We label the graphs
with $\mathcal{R}_{M}$ and $\mathcal{R}_{J}$ in a manner analogous to Eq. (56)
and Fig. 6, and use these quantities to measure the long-term stability of the
code. Their measured values are consistent with the ones for the central
density ($\mathcal{R}$). The behaviour of $\mathcal{R}_{M}$ is quite similar
for the 3 different grid resolutions, but $R_{J}$ shows considerable
dependence on the grid resolution. We have seen above (see Figure 6) that the
central density fluctuation is significantly dependent on grid resolution only
for the rotating models. We conclude that the main reason for this resolution-
sensitive behavior is the fact that angular momentum conservation is sensitive
to grid resolution. Therefore, simulations for rotating stars require high
grid resolution, otherwise angular momentum conservation will fail, and other
stationary properties of the start (such as central density) will also show
substantial, and non-physical, time evolution.

##  6 Radial Pulsation Frequency Test

Even without any explicitly-added perturbations, it is natural for our
numerical simulation of stationary stars to give rise to normal mode
oscillations due to intrinsic numerical errors. These errors occur for a
variety of reasons, including 1) truncation error due to the discretization
scheme, 2) the artificial atmosphere (floor) whereby the primitive variables
(pressure, density) are restricted from falling below minimum values to avoid
code crashes (the sound velocity becomes unbounded when vacuum is encountered
in the numerical calculations), and 3) the numerical limitation in describing
the stellar surface. Furthermore, the artificial atmosphere is known to excite
higher overtone modes.

The frequencies of various modes depend only on the structure of a given star,
and can be calculated by various methods. As explained above, our stationary
models oscillate even when we do not explicitly introduce external or internal
perturbations. We attempted to compare the frequencies of the modes excited in
our models with those obtained by normal mode analysis. The fundamental mode
(F-mode hereafter) frequency is very closely related to the dynamical time
($\sim 1/\sqrt{\rho}$) and the associated overtones have frequencies of
similar order.

Although using calculations based on cylindrical coordinates is not an
efficient way to compute radial pulsations, our code should still be able to
approximately compute the correct pulsation frequencies. The detailed
perturbation formulations and numerical methods we use for investigating the
radial pulsations are described in Appendix A. For initial conditions we use a
non-rotating equilibrium star with a baryon mass $1.4\mbox{$\rm{M}_{\odot}$}$.
we performed the test with and without the Cowling approximation, and In order
to obtain the mode frequency from the simulations, we analyzed the fluctuation
of the maximum density with time.

Specifically, we carried out Fourier transformation on the maximum density
using the FFTW package (Frigo & Johnson, 2005). To obtain better resolution in
the frequency domain, we use the zero-padding method which adds additional
zeros at the end of the time series data, effectively using interpolation
between points following the basic Fourier transformations. During the process
of obtaining a frequency having a maximum sinusoidal amplitude, leakage may
also cause additional errors. To reduce the effects of this leakage, we
multiply the time series by a window function. Here we used the Hamming window
function defined by

$w_{j}=0.54+0.46\cos\left(\frac{2\pi j}{N}\right),$ (59)

where $j$ is the index of the grid points and $N$ is the total number of
points, prior to zero-padding (Harris, 1978).

Although, as described above, some modes are excited simply due to numerical
error, their amplitudes are too small to be accurately extracted from the
simulation. We therefore introduce an explicit perturbation which can more
strongly excite the radial modes. The perturbation that we used is

$\delta\rho_{0}=B_{s}\sin\left(\pi\frac{r}{r_{s}}\right),$ (60)

where $B_{s}$ is the perturbation amplitude which we set to $B_{s}=0.001$.

00.20.40.60.8100.20.40.60.811.21.41.6 $A/A_{\rm{F-mode}}$
$f/\sqrt{\rho_{0}^{\rm{max}}}$$F$$H_{1}$$H_{2}$$H_{3}$$H_{4}$$H_{5}$
---
00.20.40.60.8100.20.40.60.811.21.41.6 $A/A_{\rm{F-mode}}$
$f/\sqrt{\rho_{0}^{\rm{max}}}$$F$$H_{1}$$H_{2}$$H_{3}$$H_{4}$$H_{5}$
Figure 8: The mode amplitudes of maximum density as a function of frequency of
the star with a baryon mass $1.4\mbox{$\rm{M}_{\odot}$}$. The vertical red
dotted lines show the frequency of the radial pulsation modes computed using
the perturbation method. The top panel shows the result when we use the
Cowling approximation, where the gravitational potential is assumed to be
fixed. In the bottom figure, we obtain the gravitational potential every few
time steps. In the both panels the 3 curves show results obtained using 3
different grid resolutions (sky blue: $257\times 257$, dark blue: $129\times
129$, and red: $65\times 65$)

Figure 8 shows the result after Fourier transformation of the time series data
given by the differences in maximum density relative to the initial time
($\rho_{0}^{\rm{max}}(t)-\rho_{0}^{\rm{max}}(t=0)$), and using calculations at
different resolutions. For comparison purposes, the vertical red lines show
the results computed from linear analysis. The mode labeled as $F$ is the
fundamental mode, while $H_{n}$ denote the $n$-th overtone radial modes. The
results shown in the figure can be summarized as follows:

  1. 1.

The most excited mode with the perturbation given by Eq. (60) is the $F$-mode.
By changing the nature of the perturbation we could make one of the overtones
the most highly excited.

  2. 2.

At low resolution, the code cannot identify high frequency modes. The reason
for this is the lack of spatial, not temporal, resolution. The eigenfunctions
describing higher overtones have large gradients near the surface which cannot
be accurately represented in the low resolution calculations.

  3. 3.

The frequency increases when we use the Cowling approximation. This is a well-
known phenomena irrespective of whether Newtonian or general relativistic
gravitation is used. This issue is discussed in more detail in the appendix.

Table 3 shows the mode frequencies computed from linear analysis as well as
the numerical simulations. Again, the stellar model is a non-rotating
spherical star of mass 1.4 $\rm{M}_{\odot}$. The relative difference between
the linear and full numerical results is listed in the last row. Here the
numerical simulations have been carried out using the highest resolution
($257\times 257$), and we list results computed with and without the Cowling
approximation.

Table 3: Comparison of mode frequencies obtained by numerical simulation ($2\pi f$) and by linear analysis ($\sigma$)  Mode | $F$ | $H_{1}$ | $H_{2}$ | $H_{3}$
---|---|---|---|---
$f/\sqrt{\rho_{0}^{\rm{max}}}$ | $0.190$ | $0.482$ | $0.734$ | $0.974$
$\sigma/2\pi/\sqrt{\rho_{0}^{\rm{max}}}$ | $0.190$ | $0.482$ | $0.733$ | $0.974$
Error(%) | $0.000$ | $0.000$ | $0.136$ | $0.000$
$|2\pi f-\sigma|/\sigma$ |  |  |  |

The frequencies we obtained from the numerical simulation with $257\times 257$
grid resolution have relative differences from those computed from linear
analysis of at most $0.1\%$. We thus conclude that the radial mode frequencies
computed from our code agree very well with the ones calculated from linear
theory. The largest difference of 0.1% was found in the second overtone (mode
$H_{2}$), while for other modes we did not find any measurable difference.

##  7 Summary and Discussion

We have developed a new hydrodynamical code which adopts a pseudo-Newtonian
treatment of the gravitational field. This code uses the so called “Valencia
formulation” for the hydrodynamical equations. From the computational
perspective, the code is modular and includes many reconstruction schemes such
as slope limiting techniques (minmod, MC, 3rd order minmod, etc.), PPM and ENO
(WENO). In 1D shock tube tests, we assessed code accuracy relative to analytic
solutions and computed convergence rates of the errors. We found that the
minmod method gives the most diffusive results, smoothing out complex features
near discontinuities. As a result it cannot be used to accurately describe
stellar surfaces, which are characterized by stiff density changes. The MC
method gives the most promising result in the shock tube test and has second
order accuracy. It can capture discontinuities very well in the pulsation mode
test, but also yields additional non-physical effects such as the excitation
of the higher order overtones near the stellar boundary. The 3minmod and PPM
methods can provide higher order accuracy and we have found that they can also
describe the stellar surface well.

In the code we also implemented 3 different flux approximation schemes: Roe,
Marquina, and HLLE. Although the results in this paper were all computed using
the HLLE approach—which is the most dissipative of the three—we have also
found that for the simulations we have considered all produce very similar
results.

In the multigrid module for computing the gravitational potential we have
implemented both second and fourth order finite-difference discretizations.
The actual value of the gravitational potential is slightly different if we
change the order of accuracy. However, the changes of maximum density in time
show very little sensitivity to the order of approximation, and we consider
the difference between the use of the second or fourth order method to be
insignificant.

In the stationary star test which is described in section 5, we evolve
equilibrium solutions describing both non-rotating and rotating stars using
our code. Our code shows stable long-time behavior of the maximum density and
conserved quantities. Based on the rates of change in the maximum density,
total mass and total angular momentum, we estimate that our code can be used
to study evolution in excess of 3,000 dynamical times with 1% error.

In the radial mode test described in section 6, modes are obtained from the
Fourier transformation of the maximum density fluctuations. We also computed
normal modes by linear analysis (see Appendix A) and found that the mode
frequencies generated by our code agree with the results from linear analysis
almost perfectly (less than 0.14%).

This code can be applied to the following astrophysical scenarios:

  1. 1.

Phenomena associated with isolated rotating neutron stars, such as
axisymmetric pulsations. Since our approach can be applied to mildly compact
stars, it is very useful to determine the amplitudes and frequencies of the
radial and non-radial modes.

  2. 2.

Accretion disks around a neutron star or black hole. It is not sufficient to
treat a disk around a compact object using Newtonian gravity, since the
gravitational field is not weak there. In addition, because the rotational
velocity of the disk is a significant fraction of $c$, we should also take
into account special relativity in our treatment of the hydrodynamics. Our
code can be a very good tool for accretion disk studies.

This work was supported by the NRF grant 2006-341-C00018, by NSERC, and by the
CIFAR Cosmology and Gravity Program.

## References

  * Arnowitt et al. (1962) Arnowitt, R., Deser, S. & Misner, C.W., 1962, Gravitation: An Introduction to Current Research, Wiley, New York.
  * Banyuls et al. (1997) Banyuls, F. and Font, J. A. and Ibanez, J. M. A. and Marti, J. M. A. & Miralles, J. A., 1997, ApJ, 476, 221.
  * Brandt (1977) Brandt, A., 1977, Math. Comp., 31, 333.
  * Centrella & Wilson (1984) Centrella, J.M., & Wilson, J.R., 1984, ApJS, 54, 229.
  * Colella & Woodward (1984) Colella, P. & Woodward, P. R., 1984, J. Comp. Phys., 54, 174.
  * Cox (1980) Cox, J. P., 1980, “Theory of Stellar Pulsation”, Princeton University Press.
  * Donat & Marquina (1996) Donat, R., & Marquina, A. 1996, J. Comp. Phys., 125, 42.
  * Donat (1998) Donat, R. 1998, J. Comp. Phys., 146, 58.
  * Einfeldt (1988) Einfeldt, B. 1988, SIAM Journal on Numerical Analysis, 25, 294.
  * Einfeldt et al. (1991) Einfeldt, B., Roe, P. L., Munz, C. D., & Sjogreen, B. 1991, J. Comp. Phys., 92, 273.
  * Font et al. (2000) Font, J. A., Miller, M., Suen, W.-M., & Tobias, M., 2000, Phys. Rev. D, 61, 044011.
  * Font (2008) Font, J. A., 2008, Living Rev. Relativity, 11, 7.
  * Frigo & Johnson (2005) Frigo, M. & Johnson, S. G., 2005, Proceedings of the IEEE, 93, 216.
  * Harris (1978) F. J. Harris, Proceedings of the IEEE, 1978, 66, 51.
  * Harten et al. (1983) Harten, A., Lax, P. D., & van Leer, B., 1983, SIAM Rev., 25, 35.
  * Hachisu (1986a) Hachisu, I., 1986a, ApJS, 61, 479.
  * Hachisu (1986b) Hachisu, I., 1986b, ApJS, 62, 461.
  * Harten et al. (1987) Harten, A. and Engquist, B. and Osher, S. & Chakravarthy, S. R., J. Comp. Phys., 71, 231.
  * Jiang & Shu (1996) Jiang, G.-S., & Shu, C.-W. 1996, J. Comp. Phys., 126, 202.
  * Kim et al. (2009) Kim, J., Kim, H. I., & Lee, H. M., 2009, MNRAS, 399, 229.
  * Liu et al. (1994) Liu, X.-D., Osher, S., & Chan, T. 1994, J. Comp. Phys., 115, 200.
  * Marti et al. (1991) Marti, J. M., Ibanez, J.M., & Miralles, J.A., 1991, Phys. Rev. D, 43, 3794.
  * Martí & Müller (2003) Martí, J. M & Müller, E., 2003, Living Rev. Relativity, 6, 7.
  * Misner et al. (1973) Misner, C. W., Thorne, K. S. & Wheeler, J. A., 1973, Gravitation, W. H. Freeman, San Francisco.
  * Nagakura et al. (2011) Nagakura, H., Ito, H., Kiuchi, K. & Yamada, S., 2011, ApJ, 731, 80.
  * Norman & Winkler (1986) Norman, M.L., & Winkler, K.-H.A., 1986, Why Ultrarelativistic Numerical Hydrodynamics is Difficult?, in Winkler, K.-H.A., and Norman, M.L., eds., Astrophysical Radiation Hydrodynamics, Proceedings of the NATO Advanced Research Workshop, Garching, Germany, August 2.13, 1982, NATO ASI Series C, vol. 188, pp. 449.
  * Roe (1981) Roe, P. L. 1981, J. Comp. Phys., 43, 357.
  * Roe (1985) Roe, P. L., 1985, Lect. Notes Appl. Math., 22, 163.
  * Sarbach & Tiglio (2012) Sarbach, O. & Tiglio, M., 2012, arXiv1203.6443.
  * Shibata (2003) Shibata, M., 2003, Phys. Rev. D, 67, 024033.
  * van Leer (1977) B. van Leer, 1977, J. Comp. Phys. 23, 276.
  * van Leer (1979) B. van Leer, 1979, J. Comp. Phys. 32, 101.
  * Wilson (1972) Wilson, J.R., 1972, ApJ, 173, 431.
  * Wilson & Mathew (2003) Wilson, J.R., & Mathews, G.J., 2003, Relativistic numerical hydrodynamics, Cambridge University Press, Cambridge, England.
  * Zhang & MacFadyen (2006) Zhang, W. and MacFadyen, A. I., ApJS, 164, 255.

##  Appendix A Perturbation equation

The eigenfrequencies and eigenfunctions of the radial pulsation of stars are
well-known in Newtonian hydrodynamics as well as in the general relativistic
case. However, the corresponding formulation has not been previously presented
for our pseudo-Newtonian approach. Here, we describe the linearized equations
that can be used to obtain eigenfrequencies and eigenfunctions of the normal
modes of spherical stars using this approximation, and following the general
relativistic framework described in Misner et al., 1973 (MTW hereafter).
First, to describe stellar oscillations—such as those occurring on the
surface—it is much more practical to use a Lagrangian description rather than
the Eulerian one adopted in section 2. The relation between the the Eulerian
and Lagrangian perturbation is (see, e.g., Cox 1980),

$\Delta f(t,r)=\delta f+f_{0}^{\prime}\zeta,$ (61)

where $\zeta$ is a Lagrangian variation in space. The law of baryon number
conservation($\nabla_{\mu}(nu^{\mu})=0$) gives

$\Delta
n=-n_{0}[r^{-2}\alpha_{0}^{3}(r^{2}\alpha_{0}^{-3}\zeta)^{\prime}-3\alpha_{0}^{-1}\delta\alpha],$
(62)

where $\alpha=\sqrt{1+2\Phi}$, $n$ is the baryon number density and ′ denotes
differentiation with respect to $r$ (See MTW Eq.(26.7)). The relation between
$n$ in Eq. (62) and $\rho_{0}$ is $\rho_{0}=m_{b}n$, where $m_{b}$ is baryon
mass and the subscript $0$ denotes the unperturbed state.

Another perturbation equation comes from the adiabatic equation of state which
offers a much easier way to find the pressure variation:

$\Gamma=\frac{n}{P}\frac{dP}{dn}.$ (63)

Since the Lagrangian variations commute with total differentiation (denoted by
$d$), Eq. (63) becomes

$\Gamma=\frac{n}{P}\frac{\Delta P}{\Delta n}.$ (64)

In addition, Eqs. (61), (62) and (64) give the following pressure variation
equation:

$\delta P=-\Gamma
P_{0}[r^{-2}\alpha_{0}^{3}(r^{2}\alpha_{0}^{-3}\zeta)^{\prime}-3\alpha_{0}^{-1}\delta\alpha]-\zeta
P_{0}^{\prime}.$ (65)

The energy conservation equation ($u_{\mu}\nabla_{\nu}T^{\mu\nu}$) gives

$\Delta\rho=\frac{\rho_{0}+P_{0}}{n_{0}}\Delta n.$ (66)

Note that $\rho_{0}$ is the energy density in the unperturbed state, rather
than the rest mass density used in the main text. Combining this with Eq.
(62), we obtain the equation for the energy density variation

$\delta\rho=-(\rho_{0}+P_{0})[r^{-2}\alpha_{0}^{3}(r^{2}\alpha_{0}^{-3}\zeta)^{\prime}-3\alpha_{0}^{-1}\delta\alpha]-\zeta\rho_{0}^{\prime}.$
(67)

The main difference here relative to the general relativistic case arises in
the computation of the perturbation of the gravitational potential. The
Poisson equation gives

$\frac{2}{r}(\alpha_{0}\delta\alpha)^{\prime}+(\alpha_{0}\delta\alpha)^{\prime\prime}=4\pi(\delta\rho+3\delta
P).$ (68)

Note that we should use only the Eulerian variation in this equation since
Lagrangian variation does not commute with partial differentiation. Eq.
(26.16) in MTW involves only first order differential equations—i.e. the
second order differentiations are rewritten in terms of first order ones. On
the other hand, in our case we cannot find any equations which can be used to
eliminate the second order differentiation. That means that we need to find
one more boundary condition to solve this equation.

Finally, the equation of motion of the fluid is obtained from the
4-acceleration ($a_{\mu}=u^{\nu}\nabla_{\nu}u_{\mu}$),

$(\rho_{0}+P_{0})\alpha_{0}^{-4}\ddot{\zeta}=-\delta
P^{\prime}-(\delta\rho+\delta P)\alpha_{0}^{-1}\alpha_{0}^{\prime}\\\
-(\rho_{0}+P_{0})(\alpha_{0}^{-1}\delta\alpha^{\prime}-\alpha_{0}^{-2}\alpha_{0}^{\prime}\delta\alpha).$
(69)

Under the assumption of the adiabatic nature of the oscillation, normal modes
are standing waves, and thus space and time variables can be separated as
follows:

$\zeta(r,t)=\xi(r)e^{i\sigma t}.$ (70)

Then, we can rewrite the equations using $\xi$ and $\sigma$,

$\displaystyle\delta P=-\Gamma
P_{0}[r^{-2}\alpha_{0}^{3}(r^{2}\alpha_{0}^{-3}\xi)^{\prime}-3\alpha_{0}^{-1}\delta\alpha]-\xi
P_{0}^{\prime}$ (71)
$\displaystyle\delta\rho=-(\rho_{0}+P_{0})[r^{-2}\alpha_{0}^{3}(r^{2}\alpha_{0}^{-3}\xi)^{\prime}-3\alpha_{0}^{-1}\delta\alpha]-\xi\rho_{0}^{\prime}$
(72) $\displaystyle(\rho_{0}+P_{0})\alpha_{0}^{-4}\sigma^{2}\zeta=\delta
P^{\prime}+(\delta\rho+\delta P)\alpha_{0}^{-1}\alpha_{0}^{\prime}$
$\displaystyle\qquad\qquad\qquad\qquad+(\rho_{0}+P_{0})(\alpha_{0}^{-1}\delta\alpha^{\prime}-\alpha_{0}^{-2}\alpha_{0}^{\prime}\delta\alpha)$
(73)

To solve Eqs. (71)–(73), we need to impose appropriate boundary conditions.
The first condition is that $\xi/r$ should be regular at the origin, and the
second one is that the pressure variation at the surface must vanish, i.e.,

$\displaystyle\frac{\xi}{r}={\rm finite~{}at~{}}r=0,$ (74)
$\displaystyle\Delta P(r=r_{s})=0.$ (75)

Unlike the general relativistic case, we cannot substitute $\delta\alpha$ and
$\delta\alpha^{\prime}$ in terms of other variations such as $\delta\rho$ and
$\delta P$. Therefore, we need an additional boundary condition for Eq. (68).
We use the properties of the gravitational potential to obtain extra
conditions. First, from the condition that the gravitational potential should
be regular at the center we obtain

$\delta\alpha^{\prime}=0{\rm~{}at~{}}r=0,.$ (76)

Second because the gravitational potential should fall off as $1/r$ beyond the
stellar surface, we have

$\delta\Phi^{\prime}+\frac{\delta\Phi}{r}=0.$ (77)

When we apply the above equation at the stellar boundary ($r=r_{s}$), we get

$\delta\alpha^{\prime}=-\frac{\delta\alpha^{2}-1}{2r\delta\alpha}{\rm~{}at~{}}r=r_{s}.$
(78)

Since Eqs. (71)–(73) and (68) are coupled, we use an iterative method to solve
them.

For the case of the Cowling approximation, which assumes that the
gravitational potential is fixed ($\delta\alpha=0$), the equations simplify
considerably:

$\displaystyle\delta P=-\Gamma
P_{0}[r^{-2}\alpha_{0}^{3}(r^{2}\alpha_{0}^{-3}\xi)^{\prime}]-\xi
P_{0}^{\prime}$ (79)
$\displaystyle\delta\rho=-(\rho_{0}+P_{0})[r^{-2}\alpha_{0}^{3}(r^{2}\alpha_{0}^{-3}\xi)^{\prime}]-\xi\rho_{0}^{\prime}$
(80) $\displaystyle(\rho_{0}+P_{0})\alpha_{0}^{-4}\sigma^{2}\zeta=\delta
P^{\prime}+(\delta\rho+\delta P)\alpha_{0}^{-1}\alpha_{0}^{\prime}$ (81)

If we compare the above equations with Eqs. (71)–(73), we observe that every
coefficient of $\delta\alpha$ is negative: therefore, as mentioned in the main
text, $\sigma$ increases when we apply the Cowling approximation.

We show the solution for $\xi/r$ for the $1.4\mbox{$\rm{M}_{\odot}$}$ star
with $K=100$ and $N=1$ with and without the Cowling approximation in Figure 9.
The $\sigma$ values corresponding to each mode are summarized in Table 3 which
appears in the main text.

00.20.40.60.811.21.41.61.800.20.40.60.81 $\tilde{\xi}/\tilde{\xi}_{s}$ $r/r_{s}$ |  -0.8-0.6-0.4-0.200.20.40.60.8100.20.40.60.81 $\tilde{\xi}/\tilde{\xi}_{s}$ $r/r_{s}$
---|---
-0.200.20.40.60.8100.20.40.60.81 $\tilde{\xi}/\tilde{\xi}_{s}$ $r/r_{s}$ |  -0.6-0.4-0.200.20.40.60.8100.20.40.60.81 $\tilde{\xi}/\tilde{\xi}_{s}$ $r/r_{s}$
-0.200.20.40.60.8100.20.40.60.81 $\tilde{\xi}/\tilde{\xi}_{s}$ $r/r_{s}$ |  -0.6-0.4-0.200.20.40.60.8100.20.40.60.81 $\tilde{\xi}/\tilde{\xi}_{s}$ $r/r_{s}$
Figure 9: Radial pulsation eigenfunction of a 1.4$\rm{M}_{\odot}$star. The
equation of state that we use is the polytropic one with $K=100$ and $N=1$. In
this figure, $\tilde{\xi}=\xi/r$ and $\tilde{\xi}_{s}=\tilde{\xi}(r=r_{s})$
where $r_{s}$ is the surface radius. The dashed and solid lines represent the
results with and without the Cowling approximation, respectively. Each panel
shows different modes (top-left ($F$), top-right ($H_{1}$), middle-left
($H_{2}$), middle- right ($H_{3}$), bottom-left ($H_{4}$) and bottom-right
($H_{5}$)) which have different oscillation frequencies.