text:
# Embedding Lagrangian Sink Particles in Eulerian Grids

## 1 Introduction

Simulations of gaseous collapse and accretion are ubiquitous in astrophysics, and all of them face a common difficulty: gravitational collapse leads to the formation of structures on length scales very small compared to the initial collapsing object. This leads to an enormous dynamic range that generally makes the full problem computationally infeasible. Simulations of star formation, for example, generally start with observed molecular cloud cores that are $`0.11.0`$ pc in size (Williams, Blitz, and McKee, 2000), while the stars that are the endpoints of the calculation are of order a solar radius ($`10^{11}`$ cm) in size – a dynamic range of $`10^7`$ in length.

Dynamic range is expensive for two distinct reasons. First, one requires enough resolution elements (cells for gridded codes, particles for gridless codes) to resolve both the largest and smallest structures in the problem. For an Eulerian code with no adaptivity, such as the widely used ZEUS package (Stone and Norman, 1992a, b; Stone, Mihalas, and Norman, 1992), increasing the linear resolution of a calculation by a factor $`f`$ requires increasing the number of cells in each dimension $`f`$, thereby increasing the total number of cells by a factor of $`f^N`$, where $`N`$ is the number of dimensions. Lagrangian approaches, such as SPH (Gingold and Monaghan, 1977; Lucy, 1977) and adaptive Eulerian approaches, such as AMR (Berger and Oliger, 1984; Berger and Collela, 1989; Bell et al., 1994), fare significantly better in this regard by following the mass and adding resolution elements only in regions of interest.

The second reason a large dynamic range is expensive is that the smallest and largest structures present in a problem often evolve on very disparate time scales, requiring a simulation to take an inordinate number of very small time steps. In explicit hydrodynamics codes this problem is embodied in the Courant condition (Richtmyer and Morton, 1967), which requires that the time step be less than the signal-crossing time of a resolution element. Increasing the linear resolution by a factor $`f`$ therefore generally requires multiplying the number of time steps by $`f`$ as well. Again, adaptive methods that allow different time steps for different resolution elements (Bate, Bonnell, and Price, 1995; Fisher et al., 2004) do somewhat better.

Even with the improvements in computational efficiency made possible by adaptivity, however, many interesting problems require more dynamic range than any code can handle. The time stepping constraint in particular is a significant barrier because, unlike additional cells, the additional time steps cannot easily be distributed on a parallel machine. Therefore simulators have introduced sinks: regions of a flow that accrete incoming material but that have no internal structure and therefore no requirements for high resolution in either time or space. Sinks provide a way to stop following collapse at a pre-chosen scale (hopefully) without damaging the rest of the calculation, thereby preventing the time step from grinding to zero and the number of resolution elements from running off to infinity.

Bate, Bonnell, and Price (1995) introduced the technique of sink particles in SPH codes. Sink particles are perfect absorbers that accrete other particles that approach within a certain distance. The accretion radius sets the smallest scale that the simulation can resolve. The boundary pressure of a sink particle is determined by extrapolation from the particles around them. While this technique has proven extremely useful and has been adopted in numerous places, SPH codes are not appropriate for all problems. To date there is no widely-used, well-tested, SPH code that includes magnetohydrodynamics or radiative transfer (although see Price and Monaghan (2003a, b) for a recent implementation of MHD with SPH that is still in the testing phase). In addition, SPH codes require significant artificial viscosity, which tends to cause artificial heating and smearing at shocks and interfaces. SPH is therefore of limited utility in problems where these features are important (Shapiro et al., 1996).

In contrast, Eulerian codes including MHD or radiative transfer are reasonably mature, and, using Godunov schemes, require much less artificial viscosity (Truelove et al., 1998). Eulerian codes have included sinks in the form of sink cells (Boss and Black, 1982). Similar to sink particles, sink cells allow mass to enter but not to leave, and their boundary pressures are found by extrapolation from neighboring cells. The disadvantage to this approach is that sink cells are fixed in the grid and are therefore inapplicable in cases where there are multiple accretors moving relative to one another (for example a binary system) or where one does not know in advance where an accretion center will form (for example in star formation in a turbulent medium).

In this paper we introduce a technique to embed a Lagrangian sink particle in an Eulerian code. This technique allows us to use Eulerian codes for cases where they are preferred, while retaining the flexibility of a moving sink center. While previous work has combined non-fluid particles with Eulerian hydrodynamics (Kravtsov, 2003), our approach is unique in that it allows the non-fluid particles to accrete from the gas, and therefore truly act as sinks. In § 2 we introduce the Eulerian AMR code and describe the method we use to embed Lagrangian sink particles within it. In § 3 we describe tests that we have done to evaluate the accuracy of this method. In discussing our test of a sink particle in the context of a disk, in § 3.4.2, we estimate the standard viscosity parameter $`\alpha `$ (Shakura and Sunyaev, 1973). We then consider a simple application of our technique in § 4: modeling the process of a flow changing from Bondi accretion to Bondi-Hoyle accretion as a result of an external shock. Finally, we discuss our conclusions in § 5.

## 2 Computational Methodology

### 2.1 The Eulerian Code

For the calculations presented in this paper we used our 3-D AMR code. The code includes hydrodynamics, gravity, and radiative transfer, but we shall refer only to the first two components in this paper. The hydrodynamics module solves the Euler equations for a compressible, multi-fluid system,

$`{\displaystyle \frac{\rho ^i}{t}}+\left(\rho ^i𝐯\right)`$ $`=`$ $`0`$ (1)
$`{\displaystyle \frac{}{t}}\left(\rho ^i𝐯\right)+\left(\rho ^i\mathrm{𝐯𝐯}\right)`$ $`=`$ $`{\displaystyle \underset{i}{}}P^i\rho ^i\varphi `$ (2)
$`{\displaystyle \frac{}{t}}\left(\rho ^ie^i\right)+\left[\left(\rho ^ie^i+P^i\right)𝐯\right]`$ $`=`$ $`\rho ^i𝐯\varphi ,`$ (3)
where $`\rho ^i`$ is the density of fluid $`i`$, $`𝐯`$ is the vector velocity (taken to be the same over all fluids), $`P^i`$ is the thermal pressure, and $`e^i`$ is the total non-gravitational energy per unit mass. We determine the potential $`\varphi `$ by solving the Poisson equation as described below. The code solves these equations using a conservative high-order Godunov scheme with an optimized approximate Riemann solver (Toro, 1997). The algorithm is second-order accurate in both space and time for smooth flows, and it provides robust treatment of shocks and discontinuities.

The gravitational module solves the Poisson equation

$$^2\varphi =4\pi G\underset{i}{}\rho ^i$$
(4)
on an adaptive grid hierarchy (as described below) to find the gravitational potential from a given density distribution. In each time step we compute the potential and then use it as a source term in the hydrodynamics equations as shown above (Truelove et al., 1998). The gravity module uses a multigrid iteration scheme to solve the linearized Poisson equation on each level of the adaptive hierarchy.

Each physics module operates within the AMR framework (Berger and Oliger, 1984; Berger and Collela, 1989; Bell et al., 1994). We discretize the problem domain onto a base, coarse level, denoted level 0. We dynamically create finer levels, numbered $`1,2,\mathrm{}n`$, nested within that coarse level as needed. The nesting process is recursive, so each fine level may contain even finer levels, providing no theoretical upper limit to the maximum resolution. In practice, limits of computational resources require that we select a maximum level of refinement allowed for a given calculation. The process for a time step is similarly recursive: one advances level $`0`$ through a single time step $`\mathrm{\Delta }t_0`$, then advances each subsequent level for the same amount of time. Each level has its own time step, and in general $`\mathrm{\Delta }t_{l+1}<\mathrm{\Delta }t_l`$, so after advancing level 0 we must advance level 1 through several steps of size $`\mathrm{\Delta }t_1`$, until it has advanced a total time $`\mathrm{\Delta }t_0`$ as well. At that point we apply a synchronization procedure to guarantee conservation of mass, momentum, and energy across the boundary between levels 0 and 1. However, each time we advance level 1 through time $`\mathrm{\Delta }t_1`$, we must advance level 2 through several steps of size $`\mathrm{\Delta }t_2`$, and so forth to the finest level present.

The sink particle framework we present below is not dependent on any of the details of the AMR framework, and may be applied equally well to fixed-grid codes. Therefore, we will suppress discussion of AMR details. However, we note here that we set our refinement criteria to guarantee that a sink particle’s accretion zone (see section 2.4.2) is always refined to the highest allowable level for a given calculation. When we refer to cell spacings and time steps in what follows, the cell spacings are always those of the finest AMR level. In addition, for simplicity we shall also suppress all discussion of multi-fluid issues.

### 2.2 Creation of Sink Particles

We may either introduce sink particles in the initial conditions for a calculation, or we may create them when necessary. Once we introduce a sink particle, we lose all knowledge of the flow in some region around it and assume that the gas within that region will continue to collapse beyond the scale resolved by our simulation. We therefore wish to introduce sinks only when the when there is good physical reason to believe that continuing the calculation without the sink will give inaccurate results, and that the gas in the vicinity of the sink is likely to continue collapsing past the scales resolved in our calculation.

Both of these conditions are met in cells that violate the Jeans criterion (Truelove et al., 1997, 1998):

$$\mathrm{\Delta }x<J\lambda _J=J\sqrt{\frac{\pi c_s^2}{G\rho }}.$$
(5)
Here, $`J`$ is a constant of order unity, $`\lambda _J`$ is the Jeans length in a cell of length $`\mathrm{\Delta }x`$, and $`c_s`$ and $`\rho `$ are the sound speed and density in the cell. Truelove et al. (1997) found that $`J=0.25`$ is sufficient to prevent artificial fragmentation in most of the problems they considered. We can also write this as a condition on the density for fixed $`\mathrm{\Delta }x`$,

$$\rho _J<J^2\frac{\pi c_s^2}{G\mathrm{\Delta }x^2}.$$
(6)
If the density in any cell does exceed $`\rho _J`$ with $`J=0.25`$, we create in the center of that cell a sink particle with mass

$$m_{\mathrm{sink}}=\left[\rho \rho _J\left(0.25\right)\right]\mathrm{\Delta }x^3,$$
(7)
so the density of the gas remaining in the cell is $`\rho _J\left(0.25\right)`$. We also transfer a proportional amount of momentum and energy from the gas to the sink particle.

Truelove et al. (1997) have shown that continuing a calculation that has violated the Jeans criterion will lead to artificial fragmentation. By creating sink particles in cells that violate the Jeans condition we prevent this from occurring. Also, since the density in the cell must have been lower in the time step before the sink particle was appeared, it follows from continuity that $`𝐯<0`$ in that cell. The fact that the cell violates the Jeans criterion indicates that its self-gravity has begun to become important. It is therefore likely that the gas in that cell will continue collapsing indefinitely. Thus, creating sinks in cells that violate the Jeans criterion meets the conditions that we create sinks only when necessary and only when we are confident that indefinite collapse is a valid approximation to the true behavior of the system. Note that these criteria are roughly analogous to those used by SPH codes (Bate, Bonnell, and Price, 1995; Bromm, Coppi, and Larson, 2002) to create sink particles: the particle must be in a region of converging flow that is gravitationally bound.

There is one cautionary note: for a cell to violate the Jeans criterion, its mass must be at least

$$m_{\mathrm{cell}}>\rho \mathrm{\Delta }x^3=\pi ^{3/2}J^3\frac{c_s^3}{\left(G^3\rho \right)^{1/2}}$$
(8)
where in the last step we eliminated $`\mathrm{\Delta }x`$ using equation (5). In comparison, the maximum mass possible for a stable, self-gravitating, isothermal object is the Bonnor-Ebert mass (Ebert, 1955; Bonnor, 1956), $`m_{\mathrm{BE}}=1.18c_s^3/\left(G^3\rho \right)^{1/2}`$. A Jeans-violating cell (for $`J=0.25`$) therefore has a minimum mass of

$$m_{\mathrm{cell}}>\frac{\pi ^{3/2}}{1.18}J^3m_{\mathrm{BE}}0.07m_{\mathrm{BE}}.$$
(9)
While the addition of neighboring cells will likely raise the total mass in the region above the Bonnor-Ebert mass, it is therefore still possible to create a sink particle in a region that is stable against collapse. However, a newly created sink particle has a very low accretion rate (see section 2.4.1). If the region is truly stable, gas will not continue to accrete onto the sink, the sink’s mass will remain very low. We thus have an ex post facto check on the validity of our sink particle creation method. If such a situation occurs, there is no way to do the calculation at the chosen level of resolution without violating the Jeans condition or creating a suspiciously small sink particle. The only choice is to redo the calculation at a higher resolution.

### 2.3 Merging Sink Particles

In a region of gravitational collapse, we often find that in a single time step a block of contiguous cells increases in density so that they all violate the Jeans condition and create sink particles. When this happens we wish to merge these particles, since allowing them to remain unmerged and possibly separate would risk a solution that contains resolution-dependent artificial fragmentation. In addition, when a sink particle is first created, gas usually continues to flow into the sink particle’s host cell and its neighbors. Since a newly-created sink particle’s mass and rate of gas accretion (see section 2.4.1) are very low, after each time step or two these cells will once again violate the Jeans condition. As a result, the code will create more sink particles, which ought to be merged into the already existing one. (This process generally continues until the sink particle is massive enough that its gas accretion rate prevents the density in the host cell from rising above the Jeans density.)

To deal with this phenomenon, at the end of each time step we group all the sink particles present in the calculation using a friends of friends (FOF) algorithm (Davis et al., 1985) with a linking length equal to the radius of the accretion zone (see section 2.4.2). We then merge all groups of particles that the FOF algorithm finds. We replace the merged group by a single particle at the center of mass of the group, and we add all particle quantities conservatively. Later on in the calculation, if two independently formed objects pass near one another for a short period, we may not wish to merge them. In this case we may temporarily reduce the radius of the accretion zone (which for technical reasons must be no larger than the merger radius) during the close passage and increase it again once the objects are sufficiently far apart. Since sink particles moving under gravity will spend very little time during a close passage, this will have a negligible effect on the overall accretion rate or final mass. This adjustment to the merger radius may be handled either manually or by an automated algorithm.

### 2.4 Accretion onto Sink Particles

#### 2.4.1 The Accretion Rate

Once it appears, a sink particle accretes gas from the surrounding cells. The gas in which the sink particle is embedded continues to evolve according to the Euler equations. Setting the accretion rate is therefore critical in cases where the flow onto the sink particle is subsonic, because the rate at which mass flows from Eulerian cells to the pressureless particle determines the amount of back-pressure opposing the accretion flow. The accretion formalism thus serves a function analogous to the extrapolation procedure used to find boundary pressure in SPH sink particle or Eulerian sink cell formalisms.

We characterize the relative importance of pressure versus gravity via the particle’s Bondi-Hoyle radius (Bondi, 1952),

$$r_{\mathrm{BH}}=\frac{GM}{v_{\mathrm{}}^2+c_{\mathrm{}}^2},$$
(10)
where $`M`$ is the particle’s mass and $`v_{\mathrm{}}`$ and $`c_{\mathrm{}}`$ are the velocity and sound speed of the gas far from the sink particle.

In the limit where the Bondi-Hoyle radius is much larger than a cell spacing, the choice of accretion rate and hence the pressure is irrelevant because the flow is supersonic near the sink particle. Even if the accretion rate is set too small, gas will flow into the sink particle’s host cell until its density is high enough to violate the Jeans condition. Once that happens, the code will create new sink particles that will immediately merge with the existing particle, thus setting an effective accretion rate higher than that given by the formula. In the opposite limit, $`r_{\mathrm{BH}}\mathrm{\Delta }x`$, the sink particle is simply a point mass moving through a uniform gas. Its gravity is relevant only on scales smaller than we resolve in the simulation. This is just the classical Bondi-Hoyle-Lyttleton problem, which has analytic solutions in the limits $`v_{\mathrm{}}c_{\mathrm{}}`$ and $`v_{\mathrm{}}c_{\mathrm{}}`$ (Hoyle and Lyttleton, 1939, 1940a, 1940b, 1940c; Bondi, 1952). In between these two limits, Ruffert (1994a) and Ruffert and Arnett (1994) give the approximate formula

$$\dot{M}=4\pi \rho _{\mathrm{}}G^2M^2\left[\frac{\lambda ^2c_{\mathrm{}}^2+v_{\mathrm{}}^2}{\left(c_{\mathrm{}}^2+v_{\mathrm{}}^2\right)^4}\right]^{1/2}=4\pi \rho _{\mathrm{}}r_{\mathrm{BH}}^2\left(\lambda ^2c_{\mathrm{}}^2+v_{\mathrm{}}^2\right)^{1/2}.$$
(11)
Here, $`\lambda `$ is a constant of order unity that depends on the equation of state of the gas. For an isothermal gas, $`\lambda =e^{3/2}/41.120`$, and we use that value throughout this work

For a real simulation with complex, turbulent flows, there is no obvious way to choose $`v_{\mathrm{}}`$. For symmetric accretion flows, the gas in the sink particle’s host cell is generally co-moving with the background at large distances. Since we have no better alternative, we therefore take $`v_{\mathrm{}}`$ to be the relative velocity of the sink particle and the gas in its host cell. Similarly, we use sound speed in the host cell for $`c_{\mathrm{}}`$. The choice of $`\rho _{\mathrm{}}`$ requires more discussion. We let ourselves be guided by the behavior we expect when simulating simple Bondi accretion. In that case, the density profile, which we denote by $`\alpha (x)\rho (x)/\rho _{\mathrm{}}`$, is the solution to a pair of coupled non-linear ordinary differential equations (Bondi, 1952). Here, $`xr/r_{\mathrm{BH}}`$ is the dimensionless radius. The density in the central cell should be $`\rho _{\mathrm{}}\alpha \left(\mathrm{\Delta }x/r_{\mathrm{BH}}\right)`$. We therefore set

$$\rho _{\mathrm{}}=\frac{\overline{\rho }}{\alpha \left(1.2\mathrm{\Delta }x/r_{\mathrm{BH}}\right)},$$
(12)
where $`\overline{\rho }`$ is the weighted mean density in the accretion region (see section 2.4.2). In the limit $`\mathrm{\Delta }xr_{\mathrm{BH}}`$, $`\alpha \left(1.2\mathrm{\Delta }x/r_{\mathrm{BH}}\right)1`$, so we recover the correct behavior for this case. We inserted the factor of 1.2 because we found that it gave improved results in the intermediate range $`\mathrm{\Delta }xr_{\mathrm{BH}}`$ (see section 3.2).

The accretion rate is reduced in the presence of rotation as described in section 2.4.2.

As a final note, Ruffert (1994b) gives a somewhat more complex formula than (11) that provides a slightly better fit to accretion rates found in their numerical simulations. However, even that formula is off by as much as 60% for isothermal flows with intermediate Mach numbers (Ruffert, 1996), and thus we decided against the additional complexity involved in implementing it. Fortunately, as we discuss in section 3.3, errors in the accretion formula tend to be self-correcting in an actual simulation, and thus the details of the accretion formula are not critical.

#### 2.4.2 The Accretion Zone

We wish the accretion rate to change smoothly as the sink particle moves across cell boundaries. Therefore we define an accretion zone around each sink particle. We set the accretion rate based on average properties in the accretion zone, and when the sink particle accretes mass, it does so from all cells within the accretion zone. In choosing the size of the accretion zone, there are two competing factors. Since the solution within the accretion zone is artificially affected by the accretion process, the larger the accretion zone, the larger the region in which we give up on the accuracy of the solution. However, in order to compute accurately the rate at which mass enters the accretion zone, we must have adequate resolution on its boundary. In addition, the size of the accretion zone will determine our ability to resolve anisotropies in the accretion flow. We define the accretion region as all the cells within a radius $`r_{\mathrm{acc}}`$ of the sink particle’s host cell. Based on experiments with different sizes, we adopt a value $`r_{\mathrm{acc}}=4\mathrm{\Delta }x`$ throughout this work. However, our implementation of the sink particle algorithm leaves the radius of the accretion region as a free parameter to be set at run time.

In cases where the particle’s Bondi-Hoyle radius is smaller than the accretion zone, it would be incorrect to set the accretion rate rate based on a uniform average of all cells, however. Therefore we define an accretion kernel with radius

$$r_K=\{\begin{array}{cc}\hfill \mathrm{\Delta }x/4:& r_{\mathrm{BH}}<\mathrm{\Delta }x/4\hfill \\ \hfill r_{\mathrm{BH}}:& \mathrm{\Delta }x/4r_{\mathrm{BH}}r_{\mathrm{acc}}/2\hfill \\ \hfill r_{\mathrm{acc}}/2:& r_{\mathrm{BH}}>r_{\mathrm{acc}}/2\hfill \end{array}.$$
(13)
Within the accretion zone we assign each cell a weight

$$w\mathrm{exp}\left(r^2/r_K^2\right),$$
(14)
where $`r`$ is the distance from the cell center to the sink particle. The minimum value of $`\mathrm{\Delta }x/4`$ for $`r_K`$ ensures that the accretion changes smoothly as the particle crosses cell boundaries even when $`r_{\mathrm{BH}}`$ is small. The maximum value of $`r_{\mathrm{acc}}/2`$ ensures that cells at the edge of the accretion region have little weight and thus there are no sudden changes in the accretion rate as cells enter or leave the accretion region. Once we have assigned weights to all cells in the accretion zone, we use a weighted average to set $`\overline{\rho }`$ in equation (12). We then compute the accretion rate for this time step from equation (11). Note that this procedure becomes undefined if the accretion zones of multiple sink particles overlap; for this reason we require that the sink particle merger radius always be greater than or equal to the accretion zone radius.

Thus far our algorithm has not included the effects of angular momentum, which may substantially reduce the accretion rate relative to the spherically symmetric case. Furthermore, in a rotating flow, low angular momentum gas along the polar axis accretes more easily than high angular momentum gas in the equatorial plane. Thus, it would be incorrect to use an accretion algorithm that accretes equally quickly from all cells regardless of their place in the rotating flow. We therefore choose a strategy that will both reduce the accretion rate in the presence of rotation and allow accretion to occur anisotropically from within the accretion zone.

To include rotation, we first divide the mass to be transferred to the sink particle (computed via equation 11) among the cells in the accretion zone so that each cell contributes an amount of mass proportional to its weight (computed via equation 14). We then divide each cell into $`8^3`$ identical point particles arranged in a uniform grid throughout the cell, each with $`1/8^3`$ the mass, momentum, and energy of the cell. For each point particle we compute its distance of closest approach to the sink particle if it were to travel on a purely ballistic trajectory while the sink particle moved at constant velocity. For a point particle with specific angular momentum $`j_{\mathrm{sp}}`$ and specific energy $`e_{\mathrm{sp}}`$ (kinetic plus gravitational) in the sink particle’s rest frame, this distance is

$$r_{\mathrm{min}}=\frac{GM}{2e_{\mathrm{sp}}}\left[1\sqrt{1+\frac{2j_{\mathrm{sp}}e_{\mathrm{sp}}}{\left(GM\right)^2}}\right].$$
(15)
If the point particle is not bound to the sink particle, so that $`e_{\mathrm{sp}}>0`$, we take $`r_{\mathrm{min}}=\mathrm{}`$. We do not want the sink particle to accrete material that has too much angular momentum to reach its “surface”. We therefore count up the number $`n`$ of point particles for which $`r_{\mathrm{min}}>\mathrm{\Delta }x/4`$ and reduce the amount of mass to be accreted from that cell by a factor $`n/8^3`$. (The factor of 4 in $`r_{\mathrm{min}}`$ is to ensure that the sink particle has an effective size smaller than the size of a cell; experimentation with different values from 0.1 to 0.5 produced no noticeable differences in behavior). For the host cell, since we cannot compute a meaningful specific angular momentum, we set $`n`$ equal to the maximum of the values of $`n`$ in the cells bordering the host cell if $`r_K\mathrm{\Delta }x/4`$, or $`n=0`$ (i.e. uninhibited accretion) if $`r_K<\mathrm{\Delta }x/4`$. Once this is done we subtract the appropriate amount of mass from each cell and add it to the sink particle. To ensure stability, we also set an absolute cap that no more than 25% of the mass may be removed from a cell in any single time step.

Next we must compute the amount of linear momentum the sink particle accretes. Since it represents an object far smaller than the grid size in the calculation, it should accrete negligible angular momentum and exert no torques on the gas. We therefore divide each cell’s momentum, taken in the sink particle’s rest frame, into components parallel and transverse to the radial vector connecting that cell to the sink particle. We reduce the cell’s radial momentum by a factor equal to the fraction of the cell’s mass that we have accreted, while we leave its transverse momentum unchanged. Thus, accretion preserves the radial velocity and the angular momentum of the gas. To ensure linear momentum conservation, we change the sink particle’s momentum by an amount equal and opposite to the total change in the momentum of the gas cells. Finally, we compute the new total energy of each gas cell by keeping the cell’s specific thermal energy constant while computing a new kinetic energy based on its new density and momentum.

We find that this accretion procedure conserves mass, momentum, and angular momentum to machine precision. However, for more discussion of angular momentum conservation see section 3.4.

### 2.5 Motion of Sink Particles

In each time step we update the position of each sink particle based on its current momentum, and we modify its momentum through accretion (as described in Section 2.4.1) and through gravity. For reasons of algorithmic speed, we handle the position and momentum update of sink particles in two steps. First, we change the momenta of sink particles due to their gravitational interactions with the the gas by an amount $`F_{\text{gas-part}}\mathrm{\Delta }t`$, where $`\mathrm{\Delta }t`$ is the time step and $`F_{\text{gas-part}}`$ is the gas-particle gravitational force. To ensure accuracy, we constrain the time step to require that

$$\mathrm{max}\left(v_{\text{part}}\right)\mathrm{\Delta }t<C\mathrm{\Delta }x,$$
(16)
where $`\mathrm{max}\left(v_{\text{part}}\right)`$ is the largest particle velocity in the calculation, $`C`$ is a constant of order unity (for our runs generally 0.5), and $`\mathrm{\Delta }x`$ is a cell spacing. This restriction is usually less stringent than the ordinary gas Courant condition. Particles have velocities comparable to the gas out of which they form, and the collapsing gas from which sink particles form is usually not the gas that has the highest velocity relative to the grid.

To compute the force on a particle, we use a Plummer law (Aarseth, 1963),

$$𝐅_{\text{gas-part}}=Gm\frac{\rho }{r^2+ϵ^2}\frac{𝐫}{r}𝑑V$$
(17)
where $`m`$ is the mass of the particle, $`\rho `$ is the gas density, $`𝐫`$ is the vector from the sink particle to a given cell or particle, and $`ϵ`$ is the softening length, which we leave as a parameter that may be set at runtime. In general, the softening length should be smaller than the size of the accretion region, to ensure that softening does not alter the rate at which gas crosses its boundary. Choosing a smaller softening length, however, increases the maximum velocity the gas will attain within the accretion region as it falls onto a sink particle, which will in turn decrease the time step due to the Courant condition. For our default choice of $`r_{\mathrm{acc}}=4\mathrm{\Delta }x`$, we therefore set a default value of $`ϵ=2\mathrm{\Delta }x`$. Since the number of particles is generally small, we compute gas-particle forces via a direct sum. In computing the force between a particle and the gas, we treat all cells except the particle’s host cell and its neighbors as point masses located at the cell center. To compute the gravitational force between a sink particle, its host cell, and its neighboring cells, we subdivide each cell into $`8^3`$ identical point particles, as described in section 2.4.2. We set the force between the cell and the sink particle equal to the sum of the forces between the sink particle and the $`8^3`$ point particles.

The second step is to update the positions and momenta of the particles including the effects of particle-particle interactions. We handle this step separately because particles may occasionally pass close to one another, within a few cell spacings. When this happens, a time step chosen via equation (16) may allow the particles to change their separation by a significant amount in a single update. In this case, a simple first- or second-order position and momentum update will not provide an accurate integration of the particle orbits. However, since computing gas updates is far more expensive than computing particle updates, we do not wish to reduce our overall time step and compute more gas updates. Instead, we integrate the particles forward through a time $`\mathrm{\Delta }t`$ using a Bulirsch-Stoer method with an adaptive time step and error control (Press et al., 1992). During this integration, we consider only particle-particle gravitational interactions, which, as with the particle-gas interactions, we compute via a direct sum. Unlike with particle-gas interactions, we do not use a softened force law for particle-particle interactions. We have tested this method by placing two particles in an extremely eccentric orbit ($`e`$ = 0.998), where the closest approach of the particles is $`\mathrm{\Delta }x/100`$, and the gas density and temperature are chosen so that accretion is negligible. In this test we found that after several orbits the semi-major axis was conserved to $`1\%`$ and the eccentricity to $`0.1\%`$.

## 3 Tests of the Methodology

### 3.1 The Collapse of an Isothermal Sphere

We have tested this method against the analytic solution for the collapse of an isothermal sphere with a density 10% above the critical value (Shu, 1977). We consider a 1 $`M_{}`$ sphere of molecular H and He mixed in the standard cosmic abundance with a sound speed of 0.18 km/s, appropriate for $`10`$ K gas. To avoid the singular initial configuration of the Shu solution, we set our initial density and velocity profile to their analytic values after a time $`t=1.3\times 10^{12}`$ s, when the expansion wave has propagated $`2.4\times 10^{16}`$ cm, or 8 cells, from the origin. There are 64 cells in the radius of the sphere, and the sphere is motionless and centered on the origin. To ensure that there are no problems when the sink particle crosses cell boundaries, we ran two more tests with identical setups except that we gave the isothermal sphere an initial uniform velocity relative to the grid. In one test we used an advection velocity of half the sound speed, and in the other twice the sound speed.

The density and velocity profiles that we find from these tests are shown in figure 1, and the mass of the sink particle versus time is shown in figure 2. The calculation reproduces the analytic solution extremely well. Errors in the velocity and density profiles are a few percent in the cells adjacent to the accretion region, dropping rapidly to $`1\%`$ as one moves further away. The advected cases show a small initial transient during which the accretion rate deviates from the theoretical value, but after a time of less than the sound-crossing time of the accretion region the flow settles into a steady state. Once this happens, the error in the accretion rate is $`1`$% in all three runs.

### 3.2 Bondi Accretion

To test how well our sink particle formalism works in cases where our resolution is marginal, so that $`r_{\mathrm{BH}}\mathrm{\Delta }x`$, we ran a series of simulations of simple Bondi accretion using a variety of values for $`r_{\mathrm{BH}}/\mathrm{\Delta }x`$. In each case we use a sphere of mixed H and He at 10 K with a radius of $`1.21\times 10^{19}`$ cm. For 32 cells across the radius of the sphere, and a 1 $`M_{}`$ central star this gives $`\mathrm{\Delta }x=r_{\mathrm{BH}}`$. We initialize the density and velocity profile of the sphere to the analytic solution for the Bondi problem and allow the calculation to run until the accretion rate reaches steady state. We then repeat the calculation with different sink particle masses and compare the accretions rates to the theoretical predictions. The results are shown in Table 1.

As the table indicates, there can be a substantial error in the accretion rate when the Bondi radius of the sink particle is comparable to a cell spacing, but the error drops off rapidly as the Bondi radius either increases or decreases. The peak error seems to occur when the Bondi radius is approximately equal to the radius of the accretion zone. For a factor of 10 difference in $`r_\mathrm{B}`$ and $`r_{\mathrm{acc}}`$ the error is $`1\%`$, while for even a factor of $`10^{1/2}`$ difference the error is only $`10\%`$ or less. We experimented with several factors of order unity in equation (12), where we set $`\rho _{\mathrm{}}`$, to see if we could decrease the error. We found that 1.2 gave the best results, but that other factors between 0.5 and 2.0 increased the error in the accretion rate by no more than $`10\%`$.

This substantial error when $`r_\mathrm{B}r_{\mathrm{acc}}`$ is not particularly surprising. The Bondi radius is the point at which the flow transitions from subsonic to supersonic. If that is equal to the accretion radius, inside which we are artificially altering the physics, then the transition is not well resolved, and the density and velocity in cells near the sink particle will have substantial errors. These errors lead the code to set an incorrect accretion rate, which in turn compounds the problem by setting a back-pressure that is too large. When the flow near the sink particle is either subsonic or supersonic, this problem does not occur and the errors are far smaller. Nonetheless, this provides an important caveat to our method: one ought not use it in cases where the flow at the sink particle is transitioning from subsonic to supersonic.

### 3.3 Bondi-Hoyle Accretion

We tested the ability of our sink particle method to handle accretion from a moving medium by simulating Bondi-Hoyle accretion. We place a 1 $`M_{}`$ sink particle in an initially uniform gas, composed of a standard interstellar mix of H and He (mean particle mass of $`2.33m_p`$) at 10 K, with a density of $`10^{25}`$ $`\text{g}\text{cm}^3`$. The gas flows past the sink particle at Mach 3. We use inflow boundary conditions in the upstream direction and outflow boundaries downstream. The resolution of the finest cells in the calculation is $`7.4\times 10^{14}\text{cm}=r_{\mathrm{BH}}/50`$, where $`r_{\mathrm{BH}}`$ is as defined by equation (10).

Using equation (11), the expected accretion rate is $`\dot{M}_{\mathrm{BH}}=1.7\times 10^{12}`$ $`M_{}\text{yr}^1`$. However, this simulation is in a regime where the interpolation formula works poorly. Ruffert (1996) performed a simulation very similar to ours (the run labeled GS, which has a small accretor, gas with polytropic index $`\gamma =1.01`$, $`=3`$), and found an accretion rate of $`2.0\times 10^{12}`$ $`M_{}\text{yr}^1=1.17\dot{M}_{\mathrm{BH}}`$, where we have scaled their dimensionless result to our dimensional parameters. Ruffert (1996) also found that both the accretion rate and the flow pattern are time-dependent, and that the flow pattern shows substantial deviations from axial symmetry on scales comparable to $`r_{\mathrm{BH}}`$. The mechanism for disrupting steady, symmetric flow is not fully understood, but Foglizzo and Ruffert (1999) suggest Rayleigh-Taylor and Kelvin-Helmholtz instabilities at the shock front as the probable cause.

Figure 3 shows the system near the end of our simulation, and Figure 4 shows the accretion rate as a function of time. As expected, the flow is time dependent and unstable, with no axial symmetry. The accretion rate requires several Bondi-Hoyle times, defined by $`t_{\mathrm{BH}}r_{\mathrm{BH}}/c_s=6.35\times 10^4`$ yr, to reach an equilibrium value, and even then it shows substantial fluctuations. The accretion rate appears to reach equilibrium after $`6t_{\mathrm{BH}}`$; the average accretion rate after that point is $`2.0\times 10^{12}M_{}\text{yr}^1=1.17\dot{M}_{\mathrm{BH}}`$, in agreement with the results of Ruffert (1996).

This test demonstrates that our method produces correct results for non-symmetric, time-dependent flows. It also illustrates an important point regarding our accretion formula, equation (11): the method is self-correcting, and thus the exact details of the accretion formula do not make much difference as long as the Bondi-Hoyle radius is well-resolved. The formula uses our best guesses for $`v_{\mathrm{}}`$ and $`\rho _{\mathrm{}}`$ based on the characteristics of the flow; however, it would be unreasonable to expect our method to correctly guess these values to the level of accuracy necessary for the formula to reproduce the Ruffert (1996) result. Instead, the accretion rate is ultimately dictated by the rate at which the ordinary, unaltered hydrodynamics of our simulation brings gas into the accretion region. If equation (11) sets an accretion rate that is too low, gas will enter the region faster than it is removed by accretion, so our guess for $`\rho _{\mathrm{}}`$ will increase and the sink particle will consume gas more quickly. The opposite effect happens if equation (11) dictates that gas be removed too quickly. Thus, our accretion rate is self-correcting.

### 3.4 Rotating Flows

#### 3.4.1 Sink Particle Results

Finally, to ensure that our method does not cause artificial accretion or angular momentum transport, we tested the evolution of a disk around a sink particle. We place a 1 $`M_{}`$ sink particle at the center of a thin 10 K gas disk with a radius $`r_0=2\times 10^{15}\text{cm}`$ and a power-law surface density profile $`\mathrm{\Sigma }=\mathrm{\Sigma }_0(r_0/r)^{k_\rho }`$, with $`\mathrm{\Sigma }_0=0.1\text{g}\text{cm}^2`$, and $`k_\rho =1`$. The disk is in Keplerian rotation, and the grid is chosen so the finest cells are $`2.1\times 10^{13}\text{cm}`$ in length. We simulated the evolution of the disk for more than 100 orbital periods at the edge of the accretion region.

In analyzing this test, it is crucial to separate the effects of the sink particle from those of ordinary numerical viscosity. As one approaches the sink particle, the circles in which the gas is flowing are resolved by fewer and fewer cells. Numerical viscosity therefore becomes significant, and will cause angular momentum transport. This effect causes evacuation of the gas in the disk near the sink particle, with some of the gas falling towards the center and the rest pushed outwards. Nelson, Benz, and Ruzmaikina (2000) and Okamoto et al. (2003) report an analogous phenomenon in SPH calculations, as do Kuznetsov et al. (1999) in a 3D Eulerian code gridded inhomogenously in spherical coordinates. To disentangle this effect from possible artificial angular momentum transport due to the sink particle, we ran two identical simulations. In one, we used the standard sink particle. In the other, we disabled accretion onto the sink particle. In this latter case, because the sink particle is far more massive than the disk and thus does not move significantly during the calculation, the sink particle’s sole effect is to impose a point gravitational potential.

Figure 5 shows the surface density versus radius at nearly identical times in the two calculations. In the simulation without accretion, there is a clear dip in the surface density between about 5 and 20 AU. Some of this gas has fallen into an unresolved hydrostatic object extending a few cells from the origin, leading to a density enhancement there. The rest has been pushed out further into the disk, leading to the alternating enhanced and diminished density between 45 and 70 AU. Examination of surface density profiles at other times shows that this phenomenon is a wave moving out from the origin, a result of material being pushed away from the sink particle by numerical viscosity.

In the calculation with sink particle accretion, the column density matches the initial surface density well except within about 20 AU of the sink particle, where it falls sharply. In this case, the gas being evacuated from the low resolution region around the sink particle has mostly been accreted. There is a slight density enhancement just outside the evacuated region, but it is smaller both in magnitude and in extent than in the calculation with no accretion. Similarly, the wave caused by material pushed outwards from the origin is far smaller in this simulation In both calculations the evacuated region is approximately the same size. The small decrease in density at radii $`>80`$ AU apparent in both simulations is a result of the pressure boundary conditions on the disk, and is unrelated to the sink particle.

Figure 6 shows the size of the evacuated region of the disk versus time. We define the edge of the evacuated region as the smallest radius for which the surface density is 90% or more of the initial surface density at that radius. In the simulation without accretion, we exclude the inner 4 cells where the hydrostatic gas has accumulated. As the plot shows, the radius of the evacuated region is slightly less in the run with sink accretion than without. This indicates that the sink particle is not causing enhanced accretion or angular momentum transport. The sink particle may actually lead to better results by preventing material that is artificially pushed outward by numerical viscosity from escaping the accretion region and affecting the rest of the calculation.

#### 3.4.2 Estimating $`\alpha `$ Due to Numerical Viscosity

As a side note, in the calculation with accretion, a power-law fit to the radius of the evacuated region versus time gives

$$\frac{r_{\mathrm{evac}}}{\mathrm{\Delta }x}=6.1\left[\frac{\mathrm{\Omega }\left(r_{\mathrm{acc}}\right)t}{2\pi }\right]^{0.23},$$
(18)
where $`\mathrm{\Omega }\left(r_{\mathrm{acc}}\right)`$ is the angular velocity of the disk at the edge of the accretion region; recall that we use $`r_{\mathrm{acc}}=4\mathrm{\Delta }x`$ in this work. We can interpret this relation as a resolution requirement for the number of cells $`r_{\mathrm{evac}}/\mathrm{\Delta }x`$ as a function of total run time $`t`$ that we need to ensure that numerical viscosity does not affect structures at a specified distance $`r_{\mathrm{evac}}`$ from the sink particle. We can also compute an approximate viscosity parameter $`\alpha `$ for this effect. For an isothermal Keplerian disk orbiting around an object of mass $`M`$, the accretion time scale at a distance $`r`$ from the central object is

$$t_{\mathrm{acc}}\frac{r^2}{\nu },$$
(19)
where $`\nu `$ is the kinematic viscosity (Lynden-Bell and Pringle, 1974). The standard $`\alpha `$ parameter is defined (Shakura and Sunyaev, 1973) so that

$$\nu =\frac{\alpha c_s^2}{\left|r\frac{d\mathrm{\Omega }}{dr}\right|}=\frac{\alpha c_s^2}{\frac{3}{2}\sqrt{\frac{GM}{r^3}}},$$
(20)
where $`c_s`$ is the sound speed in the disk and $`\mathrm{\Omega }`$ is the disk angular velocity. Eliminating $`\nu `$,

$$\alpha \frac{3}{2}\frac{\left(GMr\right)^{1/2}}{t_{\mathrm{acc}}c_s^2}=\frac{3}{2}\frac{GM}{rc_s^2}\frac{1}{\mathrm{\Omega }t_{\mathrm{acc}}}.$$
(21)
We can identify the evacuation radius and time of equation (18) with the accretion radius and timescale of equation (21). Doing so, we find that the numerical viscosity as a function of radius is

$$\alpha 78\frac{GM}{rc_s^2}\left(\frac{r}{\mathrm{\Delta }x}\right)^{2.85}=78\frac{r_\mathrm{B}}{\mathrm{\Delta }x}\left(\frac{r}{\mathrm{\Delta }x}\right)^{3.85},$$
(22)
where $`r_\mathrm{B}GM/c_s^2`$ is the standard Bondi radius. Note that the exponent of the relation for $`\alpha `$ is $`2.85`$ rather than $`4.35=1/0.23`$ because $`\mathrm{\Omega }r^{3/2}`$. To give some feel for the consequences of this relation, for $`r_\mathrm{B}/\mathrm{\Delta }x=10`$ the numerical viscosity drops to $`\alpha =0.01`$ at $`r=18.6\mathrm{\Delta }x`$.

Equation (22) provides a useful resolution requirement for hydrodynamics codes involving disks, as one can only believe the results of a disk simulation in those regions with resolution high enough that the effective $`\alpha `$ due to numerical viscosity is much smaller than the $`\alpha `$ that arises from whatever sources of physical viscosity are present. The scaling of $`\alpha `$ with $`r`$, $`r_\mathrm{B}`$, and $`\mathrm{\Delta }x`$ depends only on geometry and on the physics of $`\alpha `$-disks, and so is likely to be about the same in any code using Cartesian geometry; the constant of proportionality, $`78`$ for our code, likely depends on the hydrodynamic algorithm. Hydrodynamic codes using different grid geometries, on the other hand, would likely have very different exponent in (22). The rather large value of $`\alpha `$ we find is for a Cartesian grid. A disk centered on the origin of a cylindrical or spherical grid would probably show much less numerical viscosity.

## 4 From Bondi to Bondi-Hoyle Accretion

### 4.1 Background

While there have been extensive studies of Bondi-Hoyle accretion for flow with uniform velocities or smooth velocity gradients (Ruffert, 1994a; Ruffert and Arnett, 1994; Ruffert, 1994b; Ruffert and Anzer, 1995; Ruffert, 1995, 1996; Foglizzo and Ruffert, 1997, 1999; Foglizzo, 2002), in many cases accretion occurs in a supersonically turbulent medium in which there are numerous shocks present. An example of accretion in a shock-filled medium is the star formation process. Observations of star forming regions show they they are turbulent, with complex morphologies and velocity structures. Non-thermal linewidths within star-forming regions range from transonic for the length scale of a core that forms a single star or small multiple system to highly supersonic on the scale of entire molecular clouds, with a power-law relation between linewidth and size in between (Larson, 1981; Ossenkopf and MacLow, 2002). Simulations of turbulence in star forming cores show that turbulent motions are able to reproduce the molecular line emissions, aspect ratios, linewidth gradients, and linewidth-size relations (Klein et al., 2003; Padoan et al., 2003).

It is therefore interesting to consider what happens when a shock rolls over an accreting object, such as a protostar in a molecular cloud clump. Initially, the accretion rate and density and velocity profiles will look like standard Bondi accretion; after long times, they will be appropriate for Bondi-Hoyle accretion. We are interested in studying the details of the transition and the time-dependence of the accretion rate.

### 4.2 The Simulation

We simulate a Mach 3 shock impacting an accreting particle. We place a sink particle with a mass of $`M_{}`$ in a gas of H and He with a temperature of 10 K. The Bondi radius is therefore $`r_\mathrm{B}=0.12`$ pc. Figures 7 and 8 show the initial configuration of the problem. We divide the computational domain into two regions. For $`x>4r_\mathrm{B}`$, the gas initially has a density and velocity profile given by the analytic solution for Bondi accretion with $`\rho _{\mathrm{}}=10^{25}\text{ g}\text{ cm}^3`$. This density is unrealistically low for a real star-forming region; we choose it because we wish to ensure that the particle does not accrete enough to substantially change its mass over the course of the simulation. We are neglecting the self-gravity of the gas, and thus, as long as the sink particle’s mass changes negligibly, the density is not a relevant parameter of the problem and may be scaled to an arbitrary value. For $`x<4r_\mathrm{B}`$, we generate post-shock conditions appropriate for an isothermal shock of Mach number $`=3`$ moving in the $`+x`$ direction into a region with density $`\rho _{\mathrm{}}`$ that is at rest. Thus, for $`x<4r_\mathrm{B}`$, we set $`\rho =^2\rho _{\mathrm{}}`$, $`𝐯=\left(1/\right)c_s\widehat{𝐱}`$, where $`c_s`$ is the sound speed.

The domain of the computation goes from $`4r_\mathrm{B}`$ to $`4r_\mathrm{B}`$ in the $`y`$ and $`z`$ directions, and from $`6r_\mathrm{B}`$ to $`10r_\mathrm{B}`$ in the $`x`$ direction. The boundary conditions are inflow/outflow in every direction. We choose refinement criteria to guarantee that the region around the sink particle is always resolved such that, at a distance $`r`$ from the sink particle, $`r/\mathrm{\Delta }x32`$. This continues to a maximum resolution of $`\mathrm{\Delta }x=r_\mathrm{B}/512=50`$ AU. Note that, for $`=3`$, $`r_\mathrm{B}=10r_{\mathrm{BH}}`$, so the maximum resolution is 51 cells per Bondi-Hoyle radius.

Before the shock reaches the sink particle, the flow is pure Bondi accretion; for the parameters we use in the simulation, the accretion rate should be $`\dot{M}_\mathrm{B}=5.9\times 10^{11}`$ $`M_{}\text{yr}^1`$. Long after the shock passes the particle, the flow should be Bondi-Hoyle accretion with $`=3`$ and a background density of $`^2\rho _{\mathrm{}}=9\times 10^{25}\text{ g}\text{ cm}^3`$. Based on our results from Section 3.3 and those of Ruffert (1996), the mean accretion rate should then be $`\dot{M}_{\mathrm{BH}}=1.8\times 10^{11}`$ $`M_{}\text{yr}^1`$.

### 4.3 Results

Figure 9 shows a series of snapshots of the simulation. As one might expect, the configuration does not change much until the shock approaches within a distance $`r_\mathrm{B}`$ of the particle. At that point, the shock starts to bow, with the part along the axis the furthest forward as it propagates into infalling gas. Panel (a) shows this effect. Panel (b) shows that as the shock passes the particle, a dense cylinder of shocked material accumulates where flows of gas swept up in the shock converge and shock again. This is the beginning of the converging streamlines and Mach cone that are characteristic of Bondi-Hoyle accretion. In panel (c), the dense cylinder is showing the first signs of the destruction of axial symmetry. In panel (d), the asymmetry is increasing. By panel (e), within one Bondi-Hoyle radius of the sink particle there is a well-developed Mach cone, which is becoming turbulent in its interior. In panel (f), the Mach cone extends the full length of the range we plot, and appears to be fully turbulent in its interior. However, in panel (f) the density of material in the Mach cone is still inflated due to the presence of material that was part of the Bondi flow and has been swept up by the shock. The densities in the cone undergo a slow decline until reaching a steady state; however, the morphology of the cone does not change further.

Figure 10 shows the accretion rate as a function of time. As the plot shows, the accretion rate is initially flat and in good agreement with the predicted value for Bondi accretion. When the shock hits the sink particle at $`t=0`$, the accretion rate immediately increases by $`50\%`$ as the particle starts accreting high density shocked material. The accretion rate then begins to fall off until it approaches its equilibrium value. After some experimentation, we found that the overall shape curve is reasonably well-fit by an exponential of the form

$$\dot{M}=\dot{M}_{\mathrm{BH}}+\left(\dot{M}_0\dot{M}_{\mathrm{BH}}\right)\mathrm{exp}\left(\frac{t}{t_{\mathrm{trans}}}\right),$$
(23)
where $`\dot{M}_0`$ is the accretion rate immediately after the shock hits the sink particle and $`t_{\mathrm{trans}}`$ is the characteristic timescale for the accretion rate to transition from Bondi to Bondi-Hoyle. Our best-fit values for this case are $`\dot{M}_0=7.7\times 10^{11}M_{}\text{yr}^1=1.3\dot{M}_\mathrm{B}=4.3\dot{M}_{\mathrm{BH}}`$ and $`t_{\mathrm{trans}}=1.1\times 10^7\text{ yr}=1.7t_\mathrm{B}=17t_{\mathrm{BH}}`$, where $`t_\mathrm{B}=r_\mathrm{B}/c_s=6.35\times 10^3`$ yr and $`t_{\mathrm{BH}}=r_{\mathrm{BH}}/c_s=6.35\times 10^4`$ yr.

The characteristic timescale for the transition is closer to the Bondi time than the Bondi-Hoyle time. This is not surprising in retrospect. Before the shock hits, gas out to a distance $`r_\mathrm{B}`$ from the sink particle is inflowing supersonically. In a Bondi-Hoyle flow, gas at distances $`r_\mathrm{B}r_{\mathrm{BH}}`$ does not develop supersonic inflow velocities because it does not spend enough time near the sink particle. However, even after it is shocked, the gas that was originally part of the Bondi flow retains its supersonic infall speed and is therefore likely to find its way down to the sink particle. Since the gas is coming from a distance of order $`r_\mathrm{B}`$, its characteristic timescale to reach the sink particle is of order $`t_\mathrm{B}`$. Once all the leftover gas that was part of the Bondi flow out to $`r_\mathrm{B}`$ has drained onto the accreting particle, the accretion rate drops down to what one would expect for pure Bondi-Hoyle flow. As a result of this effect, the particle accretes an additional amount of mass

$$\mathrm{\Delta }M\left(\dot{M}\dot{M}_{\mathrm{BH}}\right)𝑑t\left(\dot{M}_0\dot{M}_{\mathrm{BH}}\right)t_{\mathrm{trans}}56\dot{M}_{\mathrm{BH}}t_{\mathrm{BH}}$$
(24)
beyond what it would have accreted if the accretion rate had instantly shifted from Bondi to Bondi-Hoyle.

## 5 Summary

We have demonstrated a new method for including moving, Lagrangian sink particles in an Eulerian hydrodynamics code. Our method shows excellent agreement with analytic results for a number of test problems, even in regimes where our formula for the accretion rate is at best a guess. In the process of testing the behavior of our sink particle in the presence of rotating flows, we develop a method to parameterize the effects of numerical viscosity in disk simulations, and use it to compute $`\alpha `$ for our Cartesian AMR code. This method provides a resolution requirement for future disk simulations.

We have also solved a simple example problem using our sink particle method, and shown that it produces useful results in that case. We have discovered one limit of our method: that sink particles can produce significant errors in the accretion rate when the flow is transitioning from subsonic to supersonic at the accretion radius. It is not clear if SPH sink particles also encounter difficulty in this regime; we have not found any tests in the literature that address the point. Regardless, with AMR one can easily avoid the regime where $`r_\mathrm{B}r_{\mathrm{acc}}`$, since once may simply increase the resolution temporarily until the sink particle accretes enough mass that $`r_\mathrm{B}>r_{\mathrm{acc}}`$.

The new technique will extend the range of Eulerian simulations, allowing them to run for longer times on problems gravitational collapse. Thus, it extends to Eulerian codes one of the heretofore unique advantages of SPH, while retaining the accurate treatment of shocks and radiative transfer capability which are found in Eulerian approaches.

The authors thank Fumitaka Nakamura and Robert Fisher for useful discussions. This work was supported by: NASA GSRP grant NGT 2-52278 (MRK); NSF grant AST-0098365 (CFM); NASA ATP grant NAG 5-12042 (CFM and RIK); the US Department of Energy at the Lawrence Livermore National Laboratory under contract W-7405-Eng-48. This research used resources of the National Energy Research Scientific Computing Center, which is supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC03-76SF00098, and the NSF San Diego Supercomputer Center through NPACI program grant UCB267.


format: text

text:
# Extinction dynamics of Lotka-Volterra ecosystems on evolving networks

## I introduction

Lotka-Volterra models of interacting species have a well established role in population ecology LV . Being inspired by an oscillatory behavior in some prey-predator systems, they are typically used to model populations on time scale shorter than lifetime of describing species. It means that long-term properties of ecosystems (macro-evolution) are usually not captured within such an approach. On the other hand, models used to describe macro-evolution very often use the dynamics that operates at the level of species rather than individuals. Such coarse-grained models usually refers to the notion of fitness of a species that is not commonly accepted BAKSNEPP .

Recently, there has been some attempts to study macro-evolution using models equipped with dynamics that operates at the level of individuals STAUFFER ; ZIA ; HALL . Taking into account that Lotka-Volterra models are relatively successful in describing many aspects of population dynamics it would be desirable to apply such an approach also to macroevolution. Some time ago Abramson introduced a discrete version of Lotka-Volterrra ecosystem ABRAMSON  and studied certain characteristics of extinctions. His model is an example of a one-dimensional food chain with $`M(100)`$ trophic levels and a single species occupying a given trophic level. Since in realistic food webs $`M46`$ with typically many species belonging to a given trophic level MCKANE ; QUINCE , these are highly nonrealistic assumptions. Nevertheless, extinction dynamics in Abramson’s model shows some features that are characteristic to Earth biosystem.

In the present paper we introduce a Lotka-Volterra model that describes a simplified ecosystem of $`N`$ species of predators and one species of preys. Our model can be thus considered as a simple food web model with only two trophic levels. Competition between predator species is described by a certain network ALBERT  of interactions whose evolution is coupled with dynamics of the model. Namely, when a certain species becomes extinct (i.e., its density falls below a certain threshold) it is replaced by new species with a newly created set of interactions with some of existing species. Despite obvious simplifications the model exhibits some properties that are typical to more complicated ecosystems, as for example power-law distributions of intervals between extinctions. Within our model we can also examine how robust this power-law distribution is. We find that under certain conditions, as for example very sparse interactions between species, or too strong dominance of a small group of species, these power-law characteristics disappear and the model is driven into a regime where extinctions have exponential distributions or where there are no extinctions and the ecosystem enters a steady state. In our opinion, such regimes might be relevant when a restricted (either in space or time) evolution of an ecosystem or its part is studied. Interestingly, a threshold value of connectivity that separates power-law extinctions and steady state is very close to the percolation threshold of the random network of inter-species interactions.

According to a large class of statistical physics models of biological evolution, avalanches of extinctions do not require external factors to trigger them, but might be a natural consequence of the dynamics of an ecosystem. As a result, these external factors, as e.g., climate changes, solar activity or impact of a big meteorite, are very often neglected in such studies ROBERTS . But such factors certainly affect the ecosystem and there is a good evidence of it NEWMAN . Let us emphasize that even the basic mechanism that triggers avalanches of extinctions is not known and is a subject of an intensive multidisciplinary debate newman\_palmer\_book .

One possibility to take external factor(s) into account in our model is to modify a growth rate of prey. Since dynamics of the model is nonlinear, such a change might have more dramatic consequences than merely a change of densities of species. And indeed we noticed that dynamics of extinctions is strongly dependent on the growth rate. It turns out, that in our model abundance of preys leads to a larger frequency of extinctions, and in periods of hunger there are less extinctions. This is clearly due to nonlinearity of the dynamics. Larger growth rate increases the density of preys that in turn increases densities of predators. With increased densities, dynamics becomes more competitive and extinctions become more frequent. Such a periodically modulated growth rate leaves some traces also in the probability distribution of extinctions. It might be interesting to notice that paleontological data also show some traces of periodic events, but their proper understanding is still missing NEWMAN ; RAUP

During evolution some species are favored and selected at the expense of less fortunate ones. Evolution constantly searches for best solutions that resembles an optimization process. For example a large size of organisms of a given species might be of advantage in some situations, but might cause some problems in the other. What might be Nature’s solution to this problem? Will it be middle-size species or rather two groups of species sitting at the extremes of conflicting requirements? In our opinion, this aspect of evolution is also often omitted in models of macroevolution. Within our model we looked at such emergent properties of species selected by evolution. It turns out that depending on some dynamical details, our model can reproduce both types of solutions of such an optimization problem.

In Section II we introduce our model and briefly describe the numerical method we used. Obtained results are presented in Section III. In Section IV we summarize our results and suggest some further extensions of our work.

## II Model and numerical calculations

We study a Lotka-Volterra ecosystem that consists of $`N`$ species of predators with densities $`\rho _i(i=1,2,\mathrm{},N)`$ who are all feeding on one species of preys with density $`\rho _0`$. We assume that each predator species $`i`$ is characterized by a parameter $`k_i`$ ($`0<k_i<1`$) that enters evolution equations of the model through death and growth terms

$`\dot{\rho _0}`$ $`=`$ $`g(t)\rho _0(1\rho _0){\displaystyle \frac{\rho _0}{N}}{\displaystyle \underset{i=1}{\overset{N}{}}}f(k_i)\rho _i`$ (1a)
$`\dot{\rho _i}`$ $`=`$ $`d(k_i)\rho _i(1\rho _0)+f(k_i)\rho _i\rho _0\left(1{\displaystyle \frac{k_i\rho _i+_j^{}k_j\rho _j}{k_i+_j^{}k_j}}\right),`$ (1b)
where $`i=1,2,\mathrm{},N`$. In our model we assume that species interact mainly through environmental capacity terms (the last term in Eq. (1b)). Namely, the growth rate of a given species $`i`$ is reduced not only due to its density but also due to weighted (with the factor $`k`$) densities of a group of randomly selected neighboring species. In Eq. (1b) summation over these neighboring species is denoted by ($`^{}`$). Approximately, we might interpret the coefficient $`k_i`$ as the size of organisms of $`i`$-th species – the bigger they are the bigger their role in the environmental capacity term. We also assume that the growth rate of preys is corrected by the environmental capacity term and due to external factors might be a slowly varying function of time ($`g(t)`$). In principle, external factors might affect also other terms of model (1), but for simplicity we restrict its influence only to the growth rate of preys. Functions $`d(k)`$ and $`f(k)`$ reflect the $`k`$-dependence of death and growth of our species. Explicit form of functions $`g(t),f(k)`$ and $`d(k)`$ will be given later.

Differential equations (1) are solved using Runge-Kutta fourth-order method. Multi-species Lotka-Volterra ecosystems were subject to intensive studies since the pioneering work of May MAY . It is known that such systems might evolve toward the steady-state with positive densities. However, in some cases, in the steady state density of some species might be zero. Each time a density of a certain species in model (1) drops below a threshold value which we fix as $`\epsilon =10^7`$ we consider such a species as extinct EPSILON . Such a species is then replaced by a new species with a randomly assigned density (from the interval (0,1)), the coefficient $`k`$ ($`0<k<1`$) that is randomly drawn from the distribution $`p(k)`$, and a new set of neighbors (all links of the ’old’ species are removed). With such rules the model rather describes $`N`$ niches, and we assume that a time to create a species that will occupy a niche is relatively short comparing to the typical lifetime of species.

We assume that a newly created species makes $`z`$ links with randomly selected neighbors. Links are not directional so a newly created species will also enter the evolution equation of species it is neighboring. If the extinct species would be chosen randomly the network of interactions would have been a random graph. However, it is the dynamics (1) that determines which species are extinct. Thus, extinct species are not selected randomly and the resulting network is in general not a random graph.

## III Results

In the following we describe numerical results obtained for some particular cases of model (1).

### III.1 Intervals between extinctions

Various paleontological data suggest that dynamics of extinctions has some power-law distributions of sizes or durations NEWMAN . In our model we measured time intervals $`t`$ between successive extinctions. In this calculations we used a constant growth term of preys $`g(t)1`$. We examined two cases: (i) model I: $`f(k_i)1,d(k_i)1`$ and (ii) model II: $`f(k_i)=k_i,d(k_i)1`$. Unless specified otherwise we select $`k_i`$ randomly with a homogeneous distribution on the interval (0,1) ($`p(k)=1`$). Our results are shown in Fig. 1. In the simplest case, model I with $`z=4`$ and $`k_i1`$ (i.e., all species during the evolution have identical $`k_i(=1)`$) we obtain exponentially decaying distribution of intervals between extinctions $`P(t)`$. Such a decay is also seen for model I (z=4) with linear distribution of $`k_i`$ namely $`p(k)=2k`$. We expect that such a behavior appears when a distribution of $`k_i`$ in the ecosystem is relatively narrow and shifted toward unity. Such an effect might be due to the small width of distribution $`p(k)`$ (i.e., a distribution from which we draw $`k_i`$) or might be dynamically generated as in model II. In this case even though $`k_i`$ are chosen from a homogeneous distribution, the dynamics favors large $`k_i`$ species (due to their larger growth rate) and they dominate the ecosystem. When the distribution of $`k_i`$ in the ecosystem is more uniform (model I with $`p(k)=1`$) our simulations suggest that $`P(t)`$ decays as a power law. Let us notice, however, that a power-law behavior is seen only on approximately one decade and we cannot exclude that on a larger time scale a different (perhaps exponential) behavior appears as was already observed in some other macroevolutionary models STAUFFER . Let us also notice that for model I with $`p(k)=k^{1/2}/2`$ the power-law distribution $`P(t)`$ seems to decay as $`t^2`$, i.e., with the exponent consistent with some paleontological data NEWMAN  as well as with predictions of some other models ZIA . However, one has to recognize that the error bars on experimental data are rather large and that a non-power law behavior cannot be excluded.

Note that a power-law decay of $`P(t)`$ is seen only for sufficiently large $`z`$. When $`z`$ is too small, we observed that the ecosystem enters the steady state where all $`\rho _i`$ are positive and there are no extinctions. This is probably due to the fact that the competition among predators is too weak (or rather too sparse). To examine the transition between this two regimes in more detail we measured the averaged time between extinctions $`\tau `$ and the results are seen in Fig. 2. One can see that $`\tau `$ diverges around $`z1.8`$ COMMENT . Such a value of the threshold parameter suggests that this transition might be related with the percolation transition in our network of interspecies interactions. To examine such a possibility we measured the average size of the largest cluster of connected links in the network $`R`$ (normalized by the number of species $`N`$) and the results are shown in Fig. 2. Vanishing of this quantity locates the percolation transition STAUFFERAHARONY . One can see that the percolation transition takes place at a larger value namely around $`z2.0`$. Our results suggest that these two transitions take place at different values of $`z`$. However the analysis of finite size effects especially in the estimation of $`\tau `$ is rather difficult and we cannot exclude that these two transitions actually overlap, as might be suggested by their proximity. Such a result would show that the dynamical regime of an ecosystem (i.e., steady state or active with power-law distribution of extinctions) is determined by the geometrical structure of its interactions.

### III.2 Effect of a modulated growth rate

Now we examine the role of a modulated in time growth rate of preys. Such a modulation is supposed to mimic the influence of an external factor like a change of a climate. One of the questions that one can ask in this context is how such a change affects the extinction dynamics. We studied model I with $`p(k)=1`$ and $`d(k_i)1`$. The growth rate of preys we chose as $`g(t)=1+A\mathrm{sin}(2\pi t/T)`$, where $`A`$ and $`T`$ are parameters. A typical behavior in case of model I with such a growth rate is shown in Fig. 3.

One can see that increased growth rate increases the density of preys $`\rho _0`$ that increases the density of predators. However, it increases also the frequency of extinctions. Such a behavior, namely increased extinction rate during abundance of food, might at first sight look as counterintuitive. This effect is related with the form of environmental capacity terms in in the growth rate in Eq. (1b), namely $`1(k_i\rho _i+_j^{}k_j\rho _j)/(k_i+_j^{}k_j)`$. Such term certainly has a larger variability for increased density of predators $`\rho _i`$, and for some species (depending on the distribution of links, coefficients $`k_i`$ and densities) it causes faster extinction. Let us also notice that since period of modulation $`T`$ is quite large, there is no retardation effect between density of preys and predators. We observed such retardation for smaller values of $`T(1000)`$.

Modulated growth rate of prays affects also the probability distribution of intervals between extinctions $`P(t)`$ as shown in Fig. 4. One can see that period of modulation $`T`$ is imprinted in $`P(t)`$. Let us notice that certain paleontological data do show some signs of periodicity but its origin still remains unclear RAUP ; NEWMAN .

It is known that slowly changing ecosystems sometimes undergo catastrophic shifts SCHEFFER . As a result, the ecosystem switches to a contrasting alternative stable state. It would be interesting to examine whether multi-species ecosystems, as described by our model (1), might also exist in such alternative states. If so, one can ask whether for example structure of the network of interspecies interactions or extinction dynamics are the same in such states.

### III.3 Emergent properties of species

It might be interesting to ask what are the characteristics of species that are preferred by the evolution in our ecosystem. Since species are characterized only by the number $`k_i`$ it is equivalent to calculating the distribution of $`k_i`$ in the steady state of model (1). Of course, due to selection this distribution in general will be different than the distribution $`p(k)`$, i.e., the distribution from which we draw $`k_i`$ of a newly created species. Some of our results are shown in Fig. 5 (all results for $`g(t)1`$, z=4, N=100).

In the case of model I ($`f(k_i)1,d(k_i)1`$) with homogeneous initial distribution of $`k_i`$ ($`p(k)=1`$) one can see that the steady state distribution is also approximately homogeneous (with a slight bias favoring small-$`k`$ species). We checked that model I shows this behavior also for other distributions $`p(k)`$ (what you put is what you get). Different behavior appears for model II ($`f(k_i)=k_i,d(k_i)1`$). In this case the growth rate factor $`f(k_i)`$ of $`i`$-th species is proportional to $`k_i`$ that certainly prefers species with large $`k_i`$. Numerical results for homogeneous distribution $`p(k)=1`$ confirm such a behavior (Fig. 5). We observed similar strong preference of large $`k_i`$ species also for Model II with other distributions $`p(k)`$.

We also examined the selection pattern in presence of some competing effects. To compensate a strong preference toward large-$`k`$ species we made simulations for our model with $`f(k_i)=k_i,d(k_i)=\sqrt{k_i}`$ and $`p(k)=1`$. Such term reduces the death rate of small-$`k`$ species. Our results show (Fig. 5) that in this case distribution of $`k_i`$ has two maxima at the extremities of the interval (0,1). On the other hand with the same model but for $`d(k_i)=(1\rho _0)^{k_i}`$ (that also reduces the death rate of small-$`k`$ species) we obtain a distribution with a single maximum around $`k=0.45`$. It would be desirable to understand the origin of the qualitative difference between these two cases.

Actually, there is yet another property of our species that is subjected to evolutionary selection, namely the number of links $`l_i`$ (degree) of a given species. Although at the beginning each species acquires $`z`$ links this number changes during the evolution because some links of a given species might be created or removed due to creation or extinction of another species. And since it is the dynamics of our model and not the random process that determines which species are removed, one can expect that the degree distribution might be different from the Poissonian distribution that is characteristic for random graphs (see BOLLO  for a precise definition of random graphs).

To check statistical properties of the network of interactions in our model we calculated the degree distribution. Our results for model I with $`z=4`$ and $`N=100`$ are shown in Fig. 6. Let us notice that although each species has $`z`$ links at the beginning it does not mean that the average number of links connected to a given site $`l_i`$ equals $`z`$ since dynamics of the model might preferentially remove sites of certain connectivity. And indeed, numerical calculations show that in this case $`l_i=2.98<z=4`$, i.e., dynamics preferentially removes sites of larger connectivity. For comparison with the random graph we plot also the Poissonian distribution $`r(l)=\mathrm{e}^{l_i}l_i^l/l!`$, where $`l_i=2.98`$. It should be emphasized that the distribution might be approximately fitted using a Poisson distribution, for example with $`l_i=2.65`$. However, it is then not a physically relevant distribution since the average connectivity $`l_i=2.65`$ differs from the value $`l_i=2.98`$ obtained from the simulations. In this sense the distribution is not Poissonian. One can see that for large connectivity the degree distribution decays faster than the Poissonian distribution. This result confirms that dynamics of the model preferentially removes highly connected species. Such sites are probably more susceptible to fluctuations in the system due to extinctions and creation of new species. On the other hand, poorly connected species are more likely to arrive at a relatively stable state. Similar results concerning the degree distribution were obtained in some other cases of our model.

## IV Conclusions

In the present paper we studied extinction dynamics of a Lotka-Volterra model of a two-level food web. In our model $`N`$ species of predators feed on a single species of preys. Competition between predators, that is specified by a certain network of interactions, leads to their extinction and replacement by new species. Distribution of intervals between successive extinctions in some cases has power-law tails and thus resembles extinction pattern of the real ecosystem. However, when the network of interactions between predators is too sparse the ecosystem enters the steady state. We have shown that such a change of behavior might be related with a percolation transition of the network. We also examined an influence of external factors on the evolution of the ecosystem. More specifically, we studied the evolution of our model in case when the growth rate of preys is changing periodically in time. It turns our that such a modulation substantially changes the frequency of extinctions. Counterintuitively, periods with abundance of preys have higher frequency of extinctions than periods with lesser amount of preys. Moreover, we examined some properties of species that are preferentially selected by the dynamics of our model. Under some conditions preferred species are a compromise to the conflicting dynamics. Under some other conditions, preferred species form two antagonistic (with respect to the conflicting rules) groups. We also examined the degree distribution of the network of interactions between species. It turns out that dynamics of the model has a slight preference to remove species of higher connectivity. As a result degree distribution shows some deviation from Poissonian distribution that is characteristic to random graphs.

It would be desirable to examine some extensions of our model. For example one can introduce additional trophic levels or other forms of interspecies interactions. One can also examine a variable number of species that would allow to create new species using certain mutation mechanism rather than assuming that they appear as soon as a niche becomes empty. Another possibility that is outside the scope of majority of macroevolutionary models would be to make further study of emergent properties of species. For example, one can imagine that a group of species in the ecosystem is well adapted and essentially not subjected to evolutionary changes. On the other hand there is a group of ’newcomers’ where evolutionary changes are much more frequent. How evolution and properties of ’newcomers’ are influenced by the properties of well-adapted species? Such problems might be easily approached within our model. Selection of a certain group of species (with a given value of $`k`$ for example) can be considered as a selection of a certain strategy. One can examine models of this kind where species are have multi-component parameters \[$`k=(k^a,k^b,\mathrm{})`$\]. Consequently, one can study evolutionary selection of more complicated traits, strategies, or behaviors. Such an approach would provide an interesting link with certain evolutionary aspects of game theory GAME .

This work was partially supported by the Swiss National Science Foundation and the project OFES 00-0578 ”COSYC OF SENS”. Some of our calculations were done on ’openMosix Cluster’ built and administrated by Lech Dȩbski at the Institute of Physics at the Adam Mickiewicz University (Poland).