Multi-scale simulations of black hole accretion in barred galaxies: Numerical methods and tests

TL;DR

This paper introduces a new numerical method for simulating gas flow in barred spiral galaxies, focusing on accurately modeling structures like nuclear rings and dust lanes using advanced computational techniques.

Contribution

It presents a novel transport scheme for angular momentum and an efficient Poisson solver for improved simulation accuracy in galactic gas dynamics.

Findings

Validated the numerical methods with tests on galactic structures

Successfully simulated nuclear rings and dust lanes

Enhanced accuracy in modeling gas flow in barred galaxies

Abstract

Due to the non-axisymmetric potential of the central bar, barred spiral galaxies form, in addition to their characteristic arms and bar, a variety of structures within the thin gas disk, like nuclear rings, inner spirals and dust-lanes. In this first of two papers, we present a method to accurately simulate the gas flow within the galactic plane in the 2D finite volume software package FOSITE, which solves the transport equations for mass, momentum and energy, and apply it to this class of objects. To this extent, we introduced a new transport scheme for angular momentum and a very efficient pseudo-spectral Poisson solver. Moreover, we provide a simple and generally applicable method of how to take care of gravity in the energy equation.

Figures

Equations

$$M\approx 1.9\cdot 10^8{\left(\frac{\sigma }{200km/s}\right)}^{5.1}{\mathrm{M}}_{\odot },$$

$$\frac{\partial \rho }{\partial t}+\nabla \cdot \left(\rho \mathbf{v}\right)=0,$$

$$\frac{\partial \left(\rho \mathbf{v}\right)}{\partial t}+\nabla \cdot \left(\rho \mathbf{v}\otimes \mathbf{v}+\mathbb{I}p\right)=0.$$

$$\frac{\partial E}{\partial t}+\nabla \cdot \left(\left(E+p\right)\mathbf{v}\right)=0$$

$$E^{\prime }=\frac{p}{\gamma -1}+\frac{1}{2}\rho {\mathbf{u}}^2.$$

$$\frac{\partial \rho }{\partial t}$$

$$\frac{\partial \left(\rho \mathbf{u}\right)}{\partial t}$$

$$\frac{\partial E^{\prime }}{\partial t}$$

$${\nabla }_R\mathbf{F}=\frac{1}{h_{\xi }{h}_{\phi }^2}\frac{\partial }{\partial \xi }\left(h_{\phi }^2{F}_1\right)+\frac{1}{h_{\xi }{h}_{\phi }}\frac{\partial }{\partial \phi }\left(h_{\xi }{F}_2\right).$$

$$\mathbf{x}=r\left(\xi \right)\left(\begin{array}{c}\cos \left(\phi \right)\\ \sin \left(\phi \right)\end{array}\right)$$

$$\frac{\partial q}{\partial t}+\nabla \cdot {\mathbf{F}}_q+\mathbf{w}\cdot \nabla q=S_q.$$

$$\frac{\partial q}{\partial t}$$

$$\mathbf{w}=r\Omega {\mathbf{e}}_{\phi }.$$

$$\mathbf{w}=\left(\tilde{w}+r\Omega \right){\mathbf{e}}_{\phi }$$

$$\frac{\mathrm{d}{u}_{i,j}}{\mathrm{d}t}=-\frac{{\mathcal{F}}_{i+\frac{1}{2}}-{\mathcal{F}}_{i-\frac{1}{2}}}{\Delta V_{i,j}}-\frac{{\mathcal{G}}_{j+\frac{1}{2}}-{\mathcal{G}}_{j-\frac{1}{2}}}{\Delta V_{i,j}}+{⟨S⟩}_{D_{i,j}}.$$

$${\mathbf{F}}_{\rho }$$

$${\mathbf{F}}_{m_{\xi }}$$

$${\mathbf{F}}_{m_{\phi }}$$

$${\mathbf{F}}_E$$

$$\frac{\mathrm{d}\rho }{\mathrm{d}t}$$

$$\frac{\mathrm{d}\left(\rho u_{\xi }\right)}{\mathrm{d}t}$$

$$\frac{\mathrm{d}\left(\rho u_{\phi }\right)}{\mathrm{d}t}$$

$$\frac{\mathrm{d}{E}^{\prime }}{\mathrm{d}t}$$

$$+\frac{w^2}{2}⟨\nabla \cdot {\mathbf{F}}_{\rho }⟩,$$

$${\nabla }^2\Phi =4\pi G\rho$$

$$\underset{|\mathbf{x}|\to \infty }{\lim }\Phi \left(\mathbf{x}\right)=0.$$

$$\Phi \left(\mathbf{x}\right)=-G{\int }_V\frac{\rho \left({\mathbf{x}}^{\prime }\right)}{\left|\mathbf{x}-{\mathbf{x}}^{\prime }\right|}d{\mathbf{x}}^{\prime }.$$

$$\frac{H\left(r\right)}{r}\ll 1.$$

$$H=H\left(r\right).$$

$$\rho \left(t,r,\phi ,z\right)=\Sigma \left(t,r,\phi \right)Z\left(r,z\right).$$

$$\Phi \left(t,r,\phi \right)$$

$$\cdot {\int }_{-\infty }^{\infty }-\frac{GZ\left(r^{\prime },z^{\prime }\right)}{\sqrt{r^2+r^{\prime 2}-2r{r}^{\prime }\cos \left(\phi -{\phi }^{\prime }\right)+z^{\prime 2}}}$$

$$\mathrm{d}{z}^{\prime }\sqrt{g}\mathrm{d}{\phi }^{\prime }\mathrm{d}{r}^{\prime }$$

$$\mathcal{G}\left(r,r^{\prime },\phi -{\phi }^{\prime }\right)={\int }_{-\infty }^{\infty }-\frac{Z\left(r^{\prime },z^{\prime }\right)}{\sqrt{r^2+r^{\prime 2}-2r{r}^{\prime }\cos \left(\phi -{\phi }^{\prime }\right)+z^{\prime 2}}}\mathrm{d}{z}^{\prime }.$$

$$Z\left(r,z\right)$$

Full text

1 Introduction

The $M$-$\sigma$-relation (Ferrarese & Merritt 2000; Gebhardt et al. 2000; McConnell et al. 2011)

[TABLE]

with $M$ the mass of the central black hole and $\sigma$ the bulge velocity dispersion, suggests a strong connection between black hole and evolution of the bulge. This correlation is unexpected, because the related length scales differ by orders of magnitudes. Black hole accretion disks release high amounts of energy and are already at comparatively moderate accretion rates very luminous (Frank et al. 2002). These accretion rates can be deduced from the luminosity of active galactic nuclei. It becomes apparent that even in the local universe a broad spectrum of accretion rates spanning $10^{-5\dots 0}\,\mathrm{M_{\odot}/yr}$ is measured (Alexander & Hickox 2012).

While the origin of the black hole-bulge connection remains generally unclear, galaxy merger (Milosavljević & Merritt 2001), AGN outflow (King 2005; King & Pounds 2015) or stellar winds (Cid Fernandes et al. 2003) seem to be connected to it. It is unknown, to which extent these processes are applicable to other galaxies and how their effect on the collective growth operates. The accretion of gas onto the black hole can only be explained by the transport of matter from much larger (kiloparsec) scales into the direct vicinity of the black hole and its gravitational sphere of influence ($<0.1\,\mathrm{pc}$). In doing so the gas not only has to loose nearly all of its angular momentum, but can also be subject to star formation by local gravitational collapse. This may suppress some of the gas flow directed into the central region. Recent high resolution observations (Xu et al. 2015; Fathi et al. 2015; Salak et al. 2016) slowly offer insight into the detailed generation of structures in the center of galaxies. Inner structures on the scale of $100\,\mathrm{pc}$ can be resolved e.g., clumpy nuclear rings, which show $1000$ times enhanced star formation rates compared to their surroundings (Xu et al. 2015) or spiral arms spread from the center, which show gas movement towards the center (Salak et al. 2016).

Often simulations of barred galaxies employ smoothed particle hydrodynamics (SPH; Shlosman 2002a; Patsis & Athanassoula 2000; Shlosman 2002b; Ann & Thakur 2005), amongst others because of the easy application and the almost exact conservation of angular momentum. A typical shortcoming of this method is the lack of resolution in the galactic center, as well as poor shock resolution, both of which, however, is important to study accretion flow onto a central black hole. Therefore in this work we use a grid-based hydrodynamic simulation code, which allows for static mesh refinement in the central region and naturally resolves shocks well. The underlying numerical scheme has been improved to ensure exact conservation of angular momentum in a rotating frame of reference, to remove any shortcomings compared to SPH.

To execute accurate measurements of the black hole accretion rate in the context of a bared galaxy, we developed a new numerical algorithm as part of the astrophysical software package FOSITE (Illenseer & Duschl 2009), which can solve systems of non-linear hyperbolic conservation equations. Hereby several computational challenges have arised, which are more generally important, not only with respect to the proposed galaxy simulations. The purpose of this paper is to introduce the new conservative transport scheme for angular momentum, the spectral self-gravity solver and the treatment of gravitational energy. The simulation of the gas flow in barred galaxies will be discussed in the second paper.

The paper is structured as follows. In section 2 we describe a method to ensure conservation of angular momentum, in section 3 a spectral self-gravity solver and in section 4 the correct treatment of gravitational energy. Section 5 is dedicated to tests of these new methods. In section 6 we conclude on our results and present an outlook on the second paper.

2 Angular momentum conservation in a rotating frame of reference

Angular momentum is one of the central conservation quantities in physics. In accretion disk physics it is of particular importance, since matter which is accreted onto the central object has to get rid if nearly all of its angular momentum. Some matter gains a considerable amount of angular momentum and migrates to outer disk parts. One of the central questions is the cause of this angular momentum and mass redistribution. To gain detailed knowledge of these processes, numerical inaccuracies have to be ruled out. Each numerical algorithm is an approximation of the analytical equations, which it is based on. Therefore small errors in the calculation are inevitable. The trick is to control the error and confine it to some well known limit. Here we want to recover the angular momentum conservation, which is guaranteed by the underlying hydrodynamical equations.

If a rotating frame of reference is taken into account, special considerations are required. A naive implementation of the resulting fictitious forces leads to a exceptionally bad angular momentum conservation (Kley 1998). Hence it requires the application of particular techniques to improve the conservation properties in this case.

Vectorial conservation laws generate, if curvilinear orthogonal coordinates are chosen, in general geometrical source terms. The cause of this is the tensor divergence (c.f. Illenseer & Duschl (2009)), which includes derivations of the basis vectors. A naive implementation would destroy the conservation properties of the vectorial conservation law.

As a generalization of a rotating reference frame general solenoidal background velocity fields are considered. In polar coordinates the transport of such velocity fields can be solved by means of a separate linear advection. Especially for quasi stationary problems the timestep limitation as implied by the Courant–Friedrichs–Lewy condition (Courant et al. 1928) of the typically supersonic rotating accretion disk can be considerably lifted, since the linear advection is unconditionally stable (Mignone et al. 2012). The original idea for this method has been proposed by Masset (2000) and is therefore often called FARGO advection. Several other astrophysical simulation codes (Mudryk & Murray 2009; Skinner & Ostriker 2010; Johnson et al. 2008; Mignone et al. 2012) have adopted the technique, which in principal accords to a moving mesh. But the mesh is only moved by multiples of complete cells and therefore the geometry remains unchanged. Since the emergence of the idea further developments of mixed Euler-Lagrangian methods have been proposed, e.g. moving annuli (Springel 2010; Duffell & MacFadyen 2012) or moving meshs of arbitrary geometry (Duffell & MacFadyen 2011). However these methods pose their own problems like higher complexity of the algorithm, a tremendous computational demand and grid noise (Duffell & MacFadyen 2012).

The common equations describing the time evolution of a compressible inviscid fluid are given by the continuity equation for the density $\rho$

[TABLE]

with the gas speed $\boldsymbol{v}$ and the momentum equation

[TABLE]

The energy equation describes the conservation of total energy, that is to say the sum of internal and kinetic energy $E=\frac{p}{\gamma-1}+\frac{1}{2}\rho\boldsymbol{v}^{2}$.

[TABLE]

If these equations are specialized to polar coordinates in a rotating reference frame, geometrical and fictitious source terms arise. Without special consideration of these in the numerical scheme, conservation of angular momentum and energy is lost.

We will now introduce a conservative implementation of the source terms. The analytical derivation is similar to the method of Paardekooper & Mellema (2006) and Mignone et al. (2012). The system can also be used in other coordinates, which are useful for rotating fluids and reference systems. Mignone et al. (2012) show derivations for spherical and cylindrical coordinates as well as the shearing box (Hawley et al. 1995).

If $\boldsymbol{v}$ is the inertial velocity field and $\boldsymbol{w}$ a solenoidal vector field, e.g. the local velocity of a rotating reference frame $\boldsymbol{w}=\boldsymbol{\Omega}\times\boldsymbol{s}$. Then $\boldsymbol{u}=\boldsymbol{v}-\boldsymbol{w}$ is called residual velocity, e.g. velocity in the rotating reference frame. Analogous the residual energy density $E^{\prime}$ is defined, which depends only on the residual velocity $\boldsymbol{u}$, but not on the inertial velocity:

[TABLE]

These definitions allow us to transform the system eq. (2) – (4) into

[TABLE]

We extend the original scheme by (Mignone et al. 2012) to arbitrary radial scalings in polar coordinates. It is therefore appropriate to use the following modified divergence operator in the azimuthal momentum equation to eliminate fictitious forces:

[TABLE]

The resulting equations and metric coefficients in polar coordinates

[TABLE]

using a monotonic differentiable radial scaling function $r\left(\xi\right)$ can be found in appendix A.

The equations Eq. (6) - Eq. (8) (and Eq. (66) - Eq. (69)) are of the type

[TABLE]

Using operator splitting (Strang 1968; McLachlan & Quispel 2002) the linear transport term $\boldsymbol{w}\cdot\nabla q$ can be removed from the rest of the equation. This yields two differential equations

[TABLE]

in which the second equation accounts for the linear transport along the direction of $\boldsymbol{w}$. In this case $\boldsymbol{w}$ points in the direction $\boldsymbol{e}_{\varphi}$ so that this equation can be solved numerically by shifting $q$ along the $\varphi$ direction. This is done by an integer translation and a simple flux transport method for the remaining part of the translation velocity $\boldsymbol{w}$. Overall an implementation of the original FARGO scheme is obtained. $\boldsymbol{w}$ is not restricted to the mean azimuthal velocity, but can accord to e.g. the velocity field of a rigid body rotation:

[TABLE]

If this is the case, the linear transport can be neglected, since it would only rotate the mesh entirely. There is no shear in the rotation motion of the mesh cell rings. Therefore this is equivalent to a rotating frame of reference with the angular velocity $\Omega$. FOSITE implements both cases so that hybrid forms

[TABLE]

are possible. The linear transport is only enabled for $\tilde{w}\neq 0$.

In case of a finite volume scheme the change of a quantity $u$ in the cell $(i,j)$ is defined by the flux of this quantity over the cell boundaries and external sources. Since FOSITE implements a semi-discrete scheme, the spatial dimensions are discretized, while the time dimension stays independent from this discretization. If we call the numerical fluxes in the different directions at the cell boundaries $\mathcal{F}$ and $\mathcal{G}$, as well as $S$ the external sources, we can write down the change of the quantity $u$ in the cell $(i,j)$ (Illenseer & Duschl 2009):

[TABLE]

In case of pure conservation laws $S=0$ holds.

Our new scheme will be conservative, if the considered quantity is exclusively modified by cell boundary fluxes. The goal is to reformulate the new source terms, which are generated by the separation of the $\boldsymbol{w}$, in terms of fluxes over the cell boundaries. The ordinary physical fluxes are identified by (Illenseer & Duschl 2009)

[TABLE]

These are subject to the classical numerical flux calculations, e.g. Riemann solvers (Harten et al. 1983; Toro et al. 1994; Roe 1981) or the Kurganov-Tadmor scheme (Kurganov & Tadmor 2000). The equations Eq. (66) - Eq. (69) are now implemented similar to the form of equation Eq. (16):

[TABLE]

in doing so terms like $\langle\nabla\cdot\boldsymbol{F}\rangle$ have to be replaced by quotient of a flux difference and the volume of a cell from equation Eq. (16). Unfortunately the derivation in Mignone et al. (2012) included minor errors, which is why we repeat it here. Note that all terms can be calculated using the original physical fluxes and the right hand side of the azimuthal momentum equation includes the complete right hand side of the continuity equation, as well as the right hand side of the energy equation includes the complete right hand side of the azimuthal momentum equation.

3 Self-gravity solver in polar coordinates

Massive accretion disks produce sufficient gravitational accelerations so that the gas mass dominates the overall gravitational potential, at least in some regions. Moreover, under certain conditions the accretion disk becomes gravitational unstable (Toomre 1964).

The Poisson equation

[TABLE]

is applied for a mass distribution with mass density $\rho$, $G$ the gravitational constant and $\Phi$ the gravitational potential. In this work we assume that $\Phi$ complies with the boundary condition of a vanishing potential at infinity

[TABLE]

Generally the Poisson equation can be solved either by a finite difference method (LeVeque 2007) or through direct integration. The first class of methods directly solves the Poisson equation Eq. (25) and requires suitable boundary conditions, which can be quite difficult.

In case of the direct integration method, we have to solve the volume integral

[TABLE]

Here the boundary conditions are not a problem, since the integral is only calculated inside the computational domain. However a difficulty arises from the evaluation at the singular point $\boldsymbol{x}^{\prime}\rightarrow\boldsymbol{x}$. Additionally, integral methods a typically computationally expensive. In the following we introduce a direct integration implemented with Fourier transformations, which can partially be cached and reused. Singular points can be avoided through introduction of a staggered mesh for the Fourier space evaluations. The method is based on Chan et al. (2006) and Li et al. (2009), but differs in the choice of the secondary grid and is generalized to polar coordinates with an arbitrary radial scaling.

We assume for simulations of accretion disks in 2D polar coordinates that they are geometrically thin. Therefore, the scale height $H$ of a accretion disk fulfills

[TABLE]

Additionally we assume that the scale height depends only on the radial coordinate but not on the azimuthal coordinate $\varphi$ or the time $t$:

[TABLE]

Then using cylindrical coordinates the vertical structure can be described by a function $Z\left(r,z\right)$. The volume density $\rho$ is therefore decomposed as a product of surface density $\Sigma$ and vertical structure $Z$:

[TABLE]

Using this in Eq. (27) and solving for $z=0$ provides:

[TABLE]

Here $\sqrt{g}$ is the Jacobian determinant. The last integral is the Green’s function

[TABLE]

Since the vertical structure $Z$ is assumed to be time independent, it can be calculated once at the beginning of the simulation. All following vertical structures have known analytic solutions to the integral. The simplest case is a razor thin disk

[TABLE]

A geometrical thin, vertically isothermal, non self-gravitating accretion disk features as vertical structure a Gaussian distribution (Lynden-Bell 1969)

[TABLE]

This applies also to vertically isothermal self-gravitating accretion disks (Illenseer & Duschl 2015). Using a scale height $H\left(r\right)$ we can calculate the vertical structure and the Green’s function:

[TABLE]

using $R^{2}=\left(r^{2}+r^{\prime 2}-2rr^{\prime}\cos\left(\varphi-\varphi^{\prime}\right)\right)/H^{2}\left(r^{\prime}\right)$ and $K_{0}$ the modified Bessel function of the second kind (see addendum B). Employing the definition

[TABLE]

we can rewrite the integral Eq. (3) as:

[TABLE]

The second integral is the convolution of $\Sigma$ and $I$. If we define $\hat{\Sigma}$ and $\hat{I}$ as the Fourier transforms of $\Sigma$ and $I$ (see addendum C) and use the convolution theorem we obtain the Fourier transform of the gravitational potential

[TABLE]

This integral can easily calculated with the rectangle method. To avoid the divergence of the solution for $r\rightarrow r^{\prime}$, we use different grids for $r$ and $r^{\prime}$. $r^{\prime}$ is evaluated at the cell centers, since $\Sigma$ is also defined at these points and therefore unnecessary reconstruction is prevented. $r$ is defined at the cell boundaries. This choice alleviates the computation of the gravitational acceleration. The following steps have to be carried out for every source term evaluation:

Calculation of the Fourier transform of $\Sigma$ 2. 2.

Evaluation of the integral Eq. (38) using the rectangle method 3. 3.

Inverse Fourier transform of $\hat{\Phi}$

As preparation we have to calculate the Fourier transform of $I$ before the simulation starts. The actual algorithm at every timestep consists only of two Fourier transformations and the summation of products (see equation Eq. (38)). This method is very fast, since for the Fourier transformation very fast and efficient algorithms are available. Typically the computation of the gravitational acceleration due to self-gravity uses about one third of the wall time of a simulation. Since 1D Fourier transformations are not efficiently parallizable, we choose to only do domain decomposition along the radial direction in rings. Furthermore, this simplifies the implementation of the orbital advection scheme, which otherwise might need additional MPI communication.

In case of very high resolution simulations low frequencies in the Fourier space may be neglected for the radial summation. This does not only reduce the amount of summation operations in the radial direction, but also the amount of communication. Tests show that the truncation error for neglecting half of the modes is typically in the lower single-digit range Li et al. (2009). This scheme is also implemented in FOSITE, but it is not used in the following simulations, since at the used resolutions the communication costs are still manageable.

The gravitational acceleration can be obtained computing the gradient of the potential

[TABLE]

which is done using finite difference approximations for the acceleration at cell centers $(i,j)$ according to

[TABLE]

4 Treatment of gravity in the energy equation

Accretion disks are approximately in a balanced state of gravitational force, pressure gradient and centrifugal force (Shakura & Sunyaev 1973; Pringle 1981). The pressure gradient in the radial direction is typically small compared to the gravitational forces and can therefore be neglected. If gravitational and centrifugal forces exactly cancel each other, the system is in a state of equilibrium. Therefore, if the radial velocity is zero, the radial momentum flux vanishes. The difficulties arising from a large imbalance of internal and total energy were already noticed by Einfeldt et al. (1991). Nonetheless their method demands a transformation of the complete numerical scheme and violates its numerous different requirements and capabilities.

We propose a new method, which has the advantage to be applicable for self-gravitating potentials. The elegant alternative method described by Jiang et al. (2013) demands a complicated modification of the conservative variables in the energy equation and thus is not readily suitable for the reformulated rothalpy equation. In contrast the new method in this work enables conservative energy transport in a rotating reference frame with self-gravitation. If the common formulation of the energy equation

[TABLE]

is used for accretion disks, a huge amount of kinetic energy has to be transported in the azimuthal direction. This result of the divergence term is mostly canceled by the potential energy source term with no net change in total energy. In the most simple case the axisymmetric potential as well as the absolute value of the velocity field feature a radial gradient. The radial gradient of the gravitational potential generates a source in the energy equation, which has to be balanced by the radial velocity gradient on the left hand side. It is extremely difficult to achieve a well balanced solution, because the flux on the left hand side is subject to reconstruction and numerical flux calculation.

Similar to the well balanced formulation of the source terms to achieve exact angular momentum conservation, the potential energy is reformulated in an analytical equivalent form, in such a way as to transport the gravitational energy together with the matter. Since the gravitational and kinetic energy have a similar absolute value, but opposite signs, they cancel each other approximately. Therefore the numerical fluxes only transport the difference to the balance of gravitational and centrifugal forces. This increases the accuracy, since minor deviations do not get lost because of truncation errors (Atkinson 2008). The gravitational source term in the energy equation is expanded by

[TABLE]

Than the energy equation can be reformulated:

[TABLE]

Now the gravitational energy is a transport term and is not handled as a common source term. The implementation follows the same concept as for the balanced formulation of the angular momentum source terms. First, the numerical fluxes of eq. 4 are calculated, than the gravitational energy is added and the spatial derivative is computed. The second term on the right hand side is handled by multiplication of the continuity equation with the potential. Thus the total energy conservation is considerably improved. All of the following simulations with energy equation and gravitational potentials are only possible because of this method. Otherwise errors are introduced typically near massive central objects in the manner of numerical instabilities, respectively negative pressure. In this case the errors in calculating internal, kinetic and potential energy are simply too large. Without the modified transport including angular momentum conservation these changes to the energy equation are from experience nonessential. Then the intrinsic dissipation of the numerical scheme is adequate to assure the stability of the solution. In the next section a simple example will test the reformulated energy equation.

5 Test problems and results

5.1 Isothermal and isentropic vortex test

The angular momentum conservation scheme is now tested with the setup of an isentropic and locally isothermal vortex. The isentropic vortex setup is based on Yee et al. (1999). The main idea is the same for both equations of state: The centrifugal forces of a rotating fluid are balanced by a radial pressure gradient. The result is a stationary flow, which has to be preserved by the numerical scheme. The original setup by Yee et al. (1999) requires to solve the energy equation. We also approximate the stationary flow in a locally isothermal setup.

Let $\rho_{\infty}$ be the background density, $T_{\infty}$ the background temperature and $u_{\infty}$, $v_{\infty}$ the background velocities in Cartesian coordinates of the undisturbed fluid. Using dimensionless quantities

[TABLE]

holds (Oertel 2005). Now we add the isentropic vortex as disturbance. An arbitrary velocity field can be chosen, which tends to zero as $r\rightarrow 0$ or $r\rightarrow\infty$. In addition the entropy shall be conserved ($\delta S=0$):

[TABLE]

Here $\beta=5$ is the vortex strength, $\gamma=1.4$ the adiabatic index and $r^{2}=x^{2}+y^{2}$. Using

[TABLE]

and the isentropic condition for an ideal gas $p/\rho^{\gamma}=p_{\infty}/\rho_{\infty}^{\gamma}=1$, we yield the analytic expression111The density solution has been corrected by a factor of $1/\pi$, which is missing in the derivation of Yee et al. (1999). for stationary flow with respect to Cartesian coordinates:

[TABLE]

If we use a locally isothermal equation of state, the density and velocity fields are defined as before. The local sound speed is a function of the density

[TABLE]

We execute simulations of the isentropic and locally isothermal vortex with and without angular momentum conservation. In addition we do these simulation in a rotating frame of reference with the angular velocity $\omega=0.79$, which corresponds to the maximum rotational velocity of the setups. At the radius of highest angular velocity the vortex has undergone about $4175$ orbits at the end of the simulation $t_{\mathrm{sim}}=2\cdot 10^{4}$. The polar grid spans a radial range from $R=0$ to $R=5$. Therefore we use AXIS conditions on the inner and NO_GRADIENTS on the outer boundary. The grid resolution is $100\times 10$ and the primitive variables are reconstructed on the cell faces using the van Leer slope limiter. The flux calculation is done with the Harten-Lax-van Leer-Contact scheme (abbr. HLLC, Toro et al. 1994).

Figure 1 shows the radial density of the isentropic, respectively locally isothermal of the simulation. The classic transport scheme is not able to preserve the vortex. After a short amount of time it has already lost most of its profile. On the contrary the modified transport scheme is able to preserve the vortex profile even for very long periods time.

Figure 2 shows the deviations from the initial configurations in the L1 norm as a function of time. The error of the simulation without modified transport is already higher than in the simulations with angular momentum conservation at any time. In the further progress the error for the classic scheme rises considerably, while it stays nearly constant for simulations with modified transport. Only at the end of the simulation some increase is noticeable.

5.2 Density-potential pairs

The first test of the new spectral self-gravity solver is the comparison of a density distribution with known analytic solution for the gravitational potential. The self-gravity module only calculates the gravitational potential of mass within the computational domain, but does no allow for boundary conditions, which account for changes of the gravitational potential due to material located outside the computational domain. Therefore density distributions should become negligible near the boundaries of the computational domain. The potential of an outside residing mass distribution can be modeled with other gravitational modules.

Useful potentials with this property are typically mass distributions, which are centered around a point in space and decline with increasing distance to this point. Such distributions are known for razor thin disks as well as disks with a Gaussian vertical structure.

In case of razor thin disks we define a around $(r_{k},\varphi_{k})$ centered mass distribution

[TABLE]

with $\sigma$ a measure for the radial extend of the disk. $R_{k}$ identifies the distance of a point in space to the center of the distribution:

[TABLE]

Such a mass distribution generates the potential (Chan et al. 2006):

[TABLE]

with $y_{k}=R_{k}/\left(2\sigma\right)$. $I_{n}$ and $K_{n}$ denote the modified Bessel functions of the first respectively second kind of integer order (see adendum B).

In the following we setup a composition of three of these mass distributions at a resolution of $1024\times 3072$ in polar coordinates, similar to a test setup of Chan et al. (2006). The solution for the gravitational potential is than the superposition of the solutions of each individual mass distribution. The surface density distribution is defined by

[TABLE]

The solution for the potential is:

[TABLE]

Figure 3 (left) shows the relative error of the numerical solution in comparison to the analytical potential. The maximum error is less than $10^{-3}$ at the points of maximum density and is typically about $10^{-5}$ if there are no large density gradients.

Next we test the implementation of the Green’s function for a vertical Gaussian distribution, using $\sigma=0.1$ as a measure for the scale height. A density distribution, which is concentrated around the $\left(r_{k},\varphi_{k}\right)$

[TABLE]

generates the gravitational potential

[TABLE]

Figure 3 (right) shows the relative error, if we use the same points in space and masses as before. The vertical Gaussian mass distribution has a maximum relative error, which is smaller than $10^{-5}$ at the points of heighest density and elsewhere of the order $10^{-7}$.

The results of density-potential-tests of Chan et al. (2006) show a similar structure of the error distribution. The maximum relative error in the simulation with vertical Gaussian mass distribution is in both works equal ( $10^{-3}$). A razor thin disk has a maximum relative error in this work of $10^{-5}$ about two orders of magnitude lower than Chan et al. (2006) results. A more detailed comparison is difficult, since Chan et al. (2006), e.g. do not document the resolution of the computational domain.

5.3 Self-gravitating rotating cylinders

To test not only the gravitational potential calculations, but also the gravitational acceleration and the behavior in a time dependent simulation, we use a setup introduced by Chan et al. (2006). In this simulation two cylinders with a Gaussian density distribution are generated, which exactly counterpart each other in the polar computational domain. A background velocity field moves the computational domain in a rigid rotation. The therefore generated centrifugal forces have to be compensated by self-gravitation of the mass distribution. This is a very demanding setup for the accuracy of the acceleration generated by the self-gravitation and the general angular momentum conservation.

The initial conditions are defined by the density distribution

[TABLE]

and the pressure

[TABLE]

as well as the velocity field

[TABLE]

Here $\operatorname{Ei}$ denotes the exponential integral (see, e.g. Abramowitz & Stegun 1965).

The simulation is executed on a polar grid with a resolution of $256\times 1024$ cells and a radial extent of $r\in\left[0.2,1.8\right]$. At these boundaries we use reflecting boundary conditions. The full set of Euler equations is closed by an isentropic equation of state with an adiabatic exponent of $\gamma=5/3$. For the flux calculation we use the Kurganov-Tadmor fluxes (Kurganov & Tadmor 2000) with linear reconstruction, slope limited by the van Leer limiter (van Leer 1974). The simulations runs for $100$ orbits.

Figure 4 shows the initial condition and the final state of the density distribution. The cylinders circular shape is very well conserved. Figure 5 shows the density in azimuthal direction at the radius of maximum density for different points in time. The density decreases slightly, but again the shape is very well conserved. The difference of the solution at the transition to the background distribution is expected, since this is an artifical cutoff.

Figure 6 shows the pressure in the first two outputs $t\in\{1,2\}$ at the beginning of the simulation. Apparently the initial condition is not in perfect equilibrium, since both cylinders eject some of their energy as a shock wave. But this is for the following examination of the conservation properties not of particular interest. Still it explains, why the decrease in peak density is stronger in the first half of the simulation.

The simulation loses in a rotating reference frame after $100$ orbits of the cylinders $0.06\%$ of the total angular momentum. Also, Chan et al. (2006) conserve the shape of the cylinders relatively well, but loose in the course of $16$ orbits already $2\%$ of the total angular momentum. Again their grid resolution is unknown, but they need $5\cdot 10^{4}$ timesteps for $16$ orbits, in contrast to our simulation, which needs less than $5\cdot 10^{5}$ for $100$ orbits. Therefore, we deduct a about $200\times$ better angular momentum conservation per timestep for our scheme with angular momentum conservation.

5.4 Vortex transport in a disk with keplerian rotation

The goal of the following test is to check the ability of the numerical scheme to preserve an anticyclonic vortex embedded in a Keplerian disk in a long-term simulation. If the preservation of a vortex is sufficiently strong, one can act on the assumption that such a vortex can be formed due to outer influences. An extensive description of the test can be found in its original source Bodo et al. (2007). Further data for comparison can be found in Mignone et al. (2012).

Initially the Keplerian background velocity is defined as $\rho_{\infty}=1$ and the pressure as $p_{\infty}=1/\left(\gamma M^{2}\right)$, using the Mach number $M=10$ at radius $r=1$. Then we add a circular vortex to the background flow:

[TABLE]

Here $x$ and $y$ are Cartesian coordinates with respect to the center of the vortex $(x_{0},y_{0})$, $\kappa=-1$ the vortex strength and $h=1/(2M)$ a dimensionless scale parameter.

We run simulations at resolutions of $256\times 1024$, $512\times 2048$ and $1024\times 4096$ with Kurganov-Tadmor (abbr. KT, Kurganov & Tadmor 2000) and HLLC fluxes. We use the same equation of state, reconstruction method and limiter as in the last section. The radial extent of the computational domain is $r\in\left[0.4,2.0\right]$. For the reflecting boundary conditions we choose a Keplerian continuation of the azimuthal velocity. The simulations run until the vortex reached $30$ orbits.

In the course of all simulations the initially circular vortex is sheared slightly along the azimuthal direction. At the same spiral density waves are generated, which slowly inch towards the domain boundaries, where they are reflected and as a consequence cross themselves. Depending on the chosen flux type and resolution, the vortex will be dissolved and the spiral arms may not be noticeable anymore.

Figure 7 shows the pressure in the simulation with HLLC fluxes without the modified gravitational energy transport scheme. A strong numerical instability emerges already after two dynamical time scales, shortly after the inner spiral arm hits the inner boundary, The instability grows quickly and results in an abort of the simulation, since the internal energy becomes negative. We observe this instability arising in simulations with angular momentum conservation independently of any other parameters or a central point mass potential, once a non-axisymmetric perturbation hits the inner boundary. Additionally its occurrence is independent from the chosen boundary condition.

Figure 8 shows the vorticity of the vortex after $30$ orbits at different resolutions. If the KT fluxes are used, the vortex is already completely destroyed at this time for the radial resolutions of $256$ and $512$. Only a radial resolution of $1024$, which is nearly inaccessible (by cpu time) in realistic, often restarted simulations, the vortex can be preserved for longer time. The simulations with HLLC fluxes show a concise anticyclonic vortex, which is to some extent already at the lowest resolution preserved. Since the start of the simulation the vortex is sheared out by the differential rotation. A inner density spiral wave can still be observed. Residuals of the outer density spiral wave, as well as its numerous reflections at the outer boundary can be adumbrated. The simulations with higher resolutions show that this one is not converged yet. A resolution of $512$ shows a similar result as the KT fluxes at $1024$. The HLLC $1024$ simulation raises the overall definition again. The vortex size diminishes and in particular the outer spiral arm is more pronounced.

Figure 9 shows the progress of the vorticity minimum for simulations of different resolutions using KT and HLLC fluxes. Figure 8 shows that after $30$ orbits the vortex is completely dissolved using KT fluxes and resolutions of $256$ and $512$. We can therefore assert that vorticies with vorticities of less than $10\%$ of their initial value, are completely dissolved. The rise in vorticity in the simulation with the coarsest resolution can be explained by accumulation of vorticity at the inner boundary, but is not related to the conservation of the vortex. The simulation with HLLC fluxes at the smallest resolution has still about $20\%$ of its initial vorticity left in the vortex center. Although the decrease in vorticity is less pronounced compared to the simulations with KT fluxes, low spatial resolution leads to a decay of the vortex likewise. This is confirmed by its change of position to smaller radii. The KT fluxes at $1024$ and the HLLC fluxes at $512$ indeed show a similar vorticity at the end of the simulation. Clearly this is the result of a switch on effect for the KT fluxes, which raises the initial vorticity at first. Afterwards the decay rate is greater in the KT fluxes simulation. The simulation with HLLC fluxes and $1024$ radial cells preserves the vortex pretty well. The vorticity raises a little bit, which can be explained by a change in the vorticity profile inside of the vortex. Of course the initial conditions do not exactly represent the by the flow adjusted solution in this setup. Overall we can note that the KT fluxes show a exponential decay at a constant rate in each simulation. In contrast the HLLC fluxes show for small resolutions a slow fading of the much weaker decay rate.

The simulations in this paper using HLLC fluxes can be compared to the FARGO simulations in Mignone et al. (2012). The flow in simulations with a resolution of $1024\times 4096$ cells show similar vortex shape and density spiral waves in both cases. The progress of the minimal vorticity is for all resolutions in very good agreement.

We can show that the vortex is dissolved within a short time, if one does not use HLLC fluxes (Toro et al. 1994) (in contrast to the HLL (Harten et al. 1983), which are very similar to the KT fluxes) and the modified formulation of the energy equation. This shows the superior definition of a HLLC solver especially for vortex flows in combination with exact angular momentum conservation.

6 Conlusions

In this paper we have introduced a new numerical algorithm to execute self-gravitating multi-scale simulations of barred galaxies with a particular interest in the black-hole accretion flow. Since angular momentum transport plays a crucial role for accretion disk dynamics, it is of great importance to reflect the angular momentum conservation in the numerical scheme. This is challenging in particular in a rotating frame of reference, since fictitious forces require careful handling. Therefore, we have introduced an implementation of a conservative numerical scheme for the solution of the Euler equations in a rotating frame of reference with angular momentum transport. This was combined with a Lagrangian advection scheme for the azimuthal transport and implemented in the astrophysical hydrodynamics software package FOSITE.

Furthermore we have implemented an efficient pseudo-spectral self-gravity solver, which is usable on polar grids with arbitrary radial scaling. The selection of a time-independent vertical disk structure allows by precalculation of a Fourier transformation to safe a considerable amount of computation time at each timestep.

At last we have introduced a modified implementation of gravitational source terms in the energy equation to avoid small but accumulating errors for the internal energy. In this method gravitational energy is advected with the flow and the gradient of the gravitational potential as source term can be avoided. The common implementation leads to non-axisymmetric instabilities arising at the inner boundary of the computational domain, especially if exact angular momentum conservation is enabled as well. We show that only the combination of angular momentum conservation, gravitational energy advection and contact resolving Riemann solvers (e.g. HLLC) is able to accurately resolve the important transport of an isentropic vortex in a Keplerian accretion disk.

In the second paper (Jung et al. 2018) we will use this new numerical algorithm to examine the self-gravitating flow of gas from galaxy regions at large radii to the smallest scales, driven by a stationary central stellar bar in a rotating frame of reference, with precise measurements of the black-hole accretion rates. This analysis will include both, isothermal simulations and simulations including the energy equation and a simple cooling mechanism to allow for marginal stable accretion disk flows.

Appendix A Euler equations in polar coordinates

The equations Eq. (6) - Eq. (8) in polar coordinates with arbitrary radial scaling are:

[TABLE]

The sources are expanded so that they can be easily implemented in dependence of the physical fluxes:

[TABLE]

The metric coefficients for polar coordinates using a monotonic differentiable radial scaling function $r\left(\xi\right)$ are:

[TABLE]

Appendix B Modified Bessel functions

The modified Bessel function of the first and second kind are (Abramowitz & Stegun 1965):

[TABLE]

Appendix C Fourier transformation

Since a number of different scalings for the Fourier transformation and convolutions are common, we note the in this work used notation and abbreviations: .

[TABLE]