Achieving Higher Order Accuracy in Space in Hydrodynamic Simulations of Self-Gravitating Gas

Modern astrophysical simulation codes employ a variety of numerical algorithms capable of achieving higher-order accuracy in both space and time. Albeit they succeed in achieving an effective higher spatial resolution and in suppressing the numerical damping of waves, to our knowledge, all current astrophysical simulations invoking self-gravity are limited to second-order accuracy in space. If we can devise an algorithm to evaluate self-gravity with a higher-order spatial accuracy, we can better the evaluation of the gravitational acceleration and gravitational energy release which dictate the evolution of many astrophysical systems. Herein, we present a numerical algorithm for self-gravitating hydrodynamics capable of achieving fourth-order accuracy for a given density distribution on a Cartesian uniform grid. First, we derive the cell-averaged gravitational potential at fourth-order accuracy from the cell-averaged density by solving the Poisson equation. Next, we obtain the cell average of the product of the density and gravitational acceleration, which differs from the cell-averaged density multiplied by the cell-averaged gravitational acceleration. We then show the verification of the algorithm by applying it to critical test problems: (1) maintaining equilibria of self-gravitating slabs, even upon advection, (2) evolving a polytropic sphere with a massive power-law envelope, and (3) conservation of specific entropy during the propagation of a sound wave.


Introduction
Gravity plays a primary role in many astronomical objects and phenomena.Thus, it is important to evaluate the gravity accurately when solving the hydrodynamic equations for modeling them.The accuracy depends on the spatial and temporal resolution and numerical scheme employed.When a numerical scheme is of n-th order accuracy in space, the numerical error is proportional to h n where h denotes the cell width in the simulation.To our knowledge, all current numerical codes available are restricted to second order accuracy in space when invoking self-gravity, even when the remainder of the physics are treated at fourth or higher order accuracy.In this, we introduce a method to evaluate the gravitational acceleration (ρg) and gravitational energy release (ρv • g) for a given density (ρ) distribution and mass flux (ρv).We show some numerical examples to demonstrate that we have achieved fouth-order accuracy in space, including (1) the evolution of a polytrope sphere and (2) the maintaining of self-gravitating slabs in equilibrium, (3) the advection of the slabs, and (4) the propagation of a sound wave modified by self-gravity.
We adopt MP5 [1] to solve the hydrodynamic equations.MP5 reconstructs cell surface values of fifth-order accuracy in space from cell averages while preserving extrema.It employs the third order Runge-Kutta method for time stepping.The temporal accuracy inherits that of the hydrodynamic equation solver employed.We use the Roe scheme as an approximate Riemann solver, although any other upwind scheme should work.
Our scheme ensures the conservation of linear momentum, i.e., the updates to momentum can be recast as the divergence of a gravitational stress tensor.As an option, we can provide a scheme to conserve total energy (including gravitational potential energy) as in Mullen et al [8].Moreover, we can evolve the the specific entropy of a gas element at third order accuracy in time, hence eliminating spurious heating or cooling exhibited by lower order schemes.
This paper is organized as follows.Section 2 describes our method and its philosophy.Our methods consist of three parts: the Poisson equation solver, the evaluation of the gravitational acceleration, and the evaluation of the gravitational energy release.Section 3 shows some numerical examples to showcase the accuracy achieved with the scheme.Section 4 summarizes the results and describes implications.

Problem Definition
We seek a fourth-order algorithm to solve the equations of self-gravitating hydrodynamics, on a uniform Cartesian grid, where Here, ρ, P , and ϕ denote the density, pressure and gravitational potential, respectively, while v and g the velocity and gravity, respectively.We assume that the flow obeys an ideal gas equation of state wherein the specific heat ratio adiabatic index γ is constant.The gravitational potential is related to the density by the Poisson equation, where G denotes the gravitational constant.
For simplicity, we ignore heating and cooling processes.Therefore, the specific entropy, s = ln P − γ ln ρ, should change in time only through numerical errors.We use this property as a proxy for the accuracy of our numerical methods.
We consider a uniform rectangular grid in Cartesian coordinates where the center of cell center is designated Here, the indices i, j, and k denote the cell number in the x-, y-, and z-directions, respectively.The grid spacing, h, is equal in all the directions.In the following, we use the half indices such as x i+1/2 = hi, y j+1/2 = hj, and z k+1/2 = hk to denote the position of the cell boundaries.

Poisson Equation
Using a fourth order accurate central difference, we discretize the Poisson Equation (Equation( 5)) as where ρ i,j,k and ϕ i,j,k denote cell averaged values.The truncation error of ϕ i,j,k is proportional to h 4 .We can solve Equation ( 7) by either multi-grid iterations [5] or the Fast Fourier Transform (FFT).Since both of these methods are established, we omit their details here.
When we apply the periodic boundary condition in all the directions, we replace the right hand side of Equation (7) with 4πG (ρ i,j,k − ρ), where ρ denotes the average density in the computation domain.We adopt this Jeans' swindle in the test problems when invoking periodic boundary conditions.

Gravitational Acceleration
Next we derive the gravitational acceleration ρg from the gravitational tensor, Since ρg = −∇ • T g , we obtain where the indices specify the components and cell surface on the right hand side.When each term on the right hand side denotes the cell surface average at fourth-order accuracy in space, the left hand side denotes the cell average gravitational acceleration of the fourth-order accuracy.Mullen et al. [8] derived the second order accurate gravitational tensor, T We can derive the other components of the tensor by changing x, y, z and indices cyclically.This second order gravitational tensor provides Since the Poisson equation is symmetric with respect to inversion, the truncation error of T (2)   g should be of the even-order.Then the fourth order accurate gravitational tensor has the form, where ∆T g denotes the correction term proportional to h 2 .We derive the correction term, ∆T g , by expanding the gravitational potential around each cell surface center up to the third order such as The cell surface averaged gravitational tensor is defined as T g,xy,i+1/2,j,k = 1 4πGh 2 y j +h/2 where we evaluate the integrand on x = x i+1/2 as a function of y and z.Substituting Equation (23) into Equations ( 24) through ( 26) we obtain 1 24 4πGT g,xz,i+1/2,j,k = ∂ϕ ∂x ∂ϕ ∂z + 1 12 We also obtain 1 24 1 24 1 24 1 24 1 24  17) through (19).Although g y,i+1,j±1/2,k , g y,i,j±1/2,k , g z,i+1,j,k±1/2 , and g z,i,j,k±1/2 have not only even-order but also odd-order truncation errors, we can eliminate the odd-order truncation errors by symmetrizing ∆T g .Thus we omit the details of the third-order truncation errors in Equations ( 31) to (34).

Gravitational Energy Release
The gravitational energy release ρv • g denotes the energy transfer rate from gravitational to other forms.The change in the gravitational energy is due to changes in the density distribution.The change in the density is described by the discretized continuity equation, where (ρv x ), (ρv y ), and (ρv z ) denote the numerical mass flux at the cell surface.Thus, it is reasonable to use the mass flux to evaluate the gravitational release.We evaluate the gravitational energy release ρv • g to be the product of the mass flux and gravity, both of which can be evaluated on the cell surface (See Equations ( 9), (10), ( 11) and ( 38)).The gravitational energy release consists of three components, The first component defined as x,i+3/2,j,k + 13 (ρv x ) i+1/2,j g x,i+1/2,j+1,k − g (4) x,i+1/2,j−1,k x,i−1/2,j+1,k − g x,i+1/2,j,k+1 − g x,i+1/2,j,k−1 x,i−1/2,j,k+1 − g (4) x,i−1/2,j,k−1 .)(40) is the fourth-order accurate cell average.We can derive the other components similarly.Note that Equation (40) takes account of the transverse gradients of mass flux and gravity in the last four lines.PPM employs a similar correction when evaluating a cell surface average (see, e.g., Equation (17) of Colella et al. [2].

Time Integration
As shown in the previous subsections, we can compute both the gravitational acceleration and gravitational energy release once the density and velocity are specified.Thus we can apply our method to each integration stage within a step when advancing the hydrodynamic equations in time.The temporal accuracy of our algorithm will then match that of the hydrodynamic scheme.Since we use a third-order Runge-Kutta integrator, the solution is of the third-order accuracy in time.

Polytrope
First we consider a polytrope in which the density and pressure are expressed as respectively, where r denotes the spherical radius from the center.Equation (41) denotes a hydrostatic equilibrium for any central density (ρ 0 ) and size scale (R).In the following numerical test, we assume ρ 0 = 3, R = 1, 4πG = 1, and γ = 5/3.We set the initial density and pressure at t = 0 to be the cell average of the analytic solution given in Equation (41).The initial velocities are set to zero.We set the mass flux on the outer boundary to zero and fix the hydrodynamic variables and the gravitational potential outside the boundary at their initial values.We compare two models having different cell widths h = 7.5 × 10 −2 and 3.75 × 10 −2 to assess the accuracy of our proposed algorithm.Figure 1 shows the relative error in the density in the z = h/2 plane of at t = 8.2.The number of numerical cells is 64 3 and 128 3 in the left and right panels, respectively.The colorbar ranges differ by a factor 16 so that the colors match when the numerical solutions are of the fourth-order accuracy.As shown in the supplementary movie, the error is due to sound waves and remains of the fourth-order accuracy during the evolution.Relatively large errors near the corner are due to the reflection of the sound waves near the boundary

Self-gravitating Sheet in Equilibrium
To quantify the numerical error, we consider the density and pressure distributions, where x = (x, y, z).The symbols k and α denote the wavenumber and amplitude.Since these profiles are periodic in all directions, we set the computational box to cover the rectangular region 0 ≤ x ≤ 2π/k x , 0 ≤ y ≤ 2π/k y , and 0 ≤ z ≤ 2π/k z .As noted in §2.2, we subtract the average density on the right hand side of Equation ( 7) since we apply periodic boundary conditions.Equations ( 42) and (43) denote a self-gravitating gas sheet in equilibrium.At the initial stage t = 0, the density and pressure are the cell average of the analytic solution.The initial velocity is zero.We use a unit system in which 4πG = 1.Any change from this state is due to numerical errors.Figure 2 shows the L1 norm of the error in the density as a function of time t for ρ 0 = 1, P 0 = 6, α = 0.3, γ = 5/3, and k = √ 6/6, √ 6/6, √ 6/3 .The three black curves denote the L1 norm for three numerical simulations computed with 32×32×16, 64×64×32, and 128×128×64 numerical cells.The red curve shows 1/16 of the L1 norm for 64 × 64 × 32 cells for comparison with that for 128 × 128 × 64 cells.The L1 norm fluctuates at early times then increases steadily at later times.The time at which the transition from fluctuations to steady increase in the error is increasingly delayed as the spatial resolution is bettered.The errors are mainly due to oscillations at early times and due to numerical diffusion in the late time.
Figure 3 is the same as Figure 2 but for the model with a initial velocity v = (0.8, 0.6, 0.0).The L1 norm fluctuates only at early times t ≤ 5.The error is larger than in the model prescribed with zero velocity.The error should be dominated by the numerical diffusion associated with advection.This model verifies not only the gravitational acceleration but also the gravitational energy release, since the velocity is large (though uniform).The same as Figure 2 but for the model with uniform velocity, v = (0.8, 0.6, 0.0).

Standing Gravito-sonic Wave
In this subsection, we consider a standing sonic wave modified by self-gravity.The initial density and pressure are assumed to be uniform at ρ 0 = 1 and P 0 = 1.We set the initial velocity to be We consider a small amplitude perturbation with α = 1.0 × 10 −3 ; we set 4πG = 1 and γ = 5/3.Equation (44) denotes a standing wave when α is infinitesimal.When α is small but finite, non-linearity of the hydrodynamics produces overtones with amplitude proportional to α 2 .The overtone modifies the density and pressure profile.However, the specific entropy can change only through numerical errors.Thus, we can measure the accuracy of the gravitational energy release by measuring the change in the specific entropy, s = ln P − γ ln ρ.
We set the initial density, momentum density (ρv), and the total energy density to be cell averaged represenations.We evaluate the specific entropy at the cell center (not the cell average) according to PPM [2].Remember that the specific entropy derived from the cell average density and pressure are not the cell average specific entropy.
Figure 4 shows the maximum specific entropy (the solid curves) and minimum (the dashed curves) for the solutions solved with 32 × 32 × 16 and 64 × 64 × 32 numerical cells.The red curve denotes 1/16 of the maximum entropy obtained with 32 × 32 × 16 cells.Figure 5 is the same as Figure 4 but for the solutions obtained with 64 × 64 × 32 and 128 × 128 × 64 cells.Since the specific entropy vanishes at the initial state, the deviation from s = 0 is of numerical origin.
Although the entropy increases with time, the absolute value is small (< 2 × 10 −8 ).Both the maximum and minimum specific entropy increase linearly as the time passes.The rate of increase is common for the maximum and minimum and proportional to cube of the cell width (ds max /dt ∝ h 3 and ds min /dt ∝ h 3 ).This numerical error is attributed to the hydrodynamic sector.The third-order convergence verifies our method to evaluate the gravitational acceleration and gravitational energy release.

Summary and Discussions
As shown in the previous sections, we have devised a numerical algorithm capable of achieving fourth-order accuracy in space in the evaluation of the gravitational acceleration ρg and the The maximum and minimum specific entropy for standing sonic waves for the solutions obtained with 32 × 32 × 16 and 64 × 64 × 32 numerical cells.The same as Figure 4 but for the solutions obtained with 64 × 64 × 32 and 128 × 128 × 64 numerical cells.gravitational energy release ρv • g.When coupling these source terms to the hydrodynamic equations, we can obtain solutions with overall fourth-order accuracy in space.The four examples shown here ensure the fourth-order accuracy in space.The fourth-order accuracy in space will improve the modeling of a self-gravitating astronomical object.
Though our scheme provides fourth-order accuracy in space, there still exist unresolved issues.First, Equations (35), (36), and (37) leave some freedoms in the defining of discretized expressions for the higher order derivatives.They might produce small (at fourth order truncation level), yet significant differences in the results (c.f., Mullen et al. [8] found that particular discretizations of the gravitational stress tensor could produce anomalous accelerations associated with ∇×g ̸ = 0).Also, we are still investigating the robustness of the algorithm against various failure modes (e.g., negative pressures); preliminary results demonstrate that we may need to reduce the spatial accuracy around shocks or introduce artificial dissipation.These areas are still under active investigation.

15thFigure 1 .
Figure 1.Relative error profiles of the density field (∆ρ/ρ), after evolving a polytrope for t ≃ 8.2 ≃ 14/ √ 4πGρ 0 The spatial resolution improves by a factor of two moving from the left to right panel.The colorbar ranges differ by a factor of 16.

15thFigure 2 .
Figure 2. The L1 norm of the error in the density as a function of time, t, for a selfgravitating gas sheet.

15th
Figure 4.The maximum and minimum specific entropy for standing sonic waves for the solutions obtained with 32 × 32 × 16 and 64 × 64 × 32 numerical cells.