A Mean-Field Approach to Simulating the Merging of Collisionless Stellar Systems Using a Particle-Based Method

TL;DR

This paper introduces a mean-field, particle-based simulation method using a self-consistent field approach to model the merging of collisionless stellar systems, demonstrating its accuracy and advantages over traditional tree codes.

Contribution

The paper develops a self-consistent field (SCF) method for simulating stellar mergers, providing a softening-free, accurate alternative to tree codes with broad application potential.

Findings

SCF and tree code results agree well on morphology and profiles.

Softening in tree codes affects orbital phase timing.

SCF method effectively models collisionless stellar mergers.

Abstract

We present a mean-field approach to simulating merging processes of two spherical collisionless stellar systems. This approach is realized with a self-consistent field (SCF) method in which the full spatial dependence of the density and potential of a system is expanded in a set of basis functions for solving Poisson's equation. In order to apply this SCF method to a merging situation where two systems are moving in space, we assign the expansion center to the center of mass of each system, the position of which is followed by a mass-less particle placed at that position initially. Merging simulations over a wide range of impact parameters are performed using both an SCF code developed here and a tree code. The results of each simulation produced by the two codes show excellent agreement in the evolving morphology of the merging systems and in the density and velocity dispersion profiles of the merged systems. However, comparing the results generated by the tree code to those obtained with the softening-free SCF code, we have found that in large impact parameter cases, a softening length of the Plummer type introduced in the tree code has an effect of advancing the orbital phase of the two systems in the merging process at late times. We demonstrate that the faster orbital phase originates from the larger convergence length to the pure Newtonian force. Other application problems suitable to the current SCF code are also discussed.

Tables

Table 1. Table 1: Initial Cartesian coordinates of the centers of the two systems.

Name	System 1	System 2
Case 0	$(-5,\mathrm{\hspace{0.17em}0},\mathrm{\hspace{0.17em}0})$	$(5,\mathrm{\hspace{0.17em}0},\mathrm{\hspace{0.17em}0})$
Case 1	$(-5,-1,\mathrm{\hspace{0.17em}0})$	$(5,\mathrm{\hspace{0.17em}1},\mathrm{\hspace{0.17em}0})$
Case 2	$(-5,-2,\mathrm{\hspace{0.17em}0})$	$(5,\mathrm{\hspace{0.17em}2},\mathrm{\hspace{0.17em}0})$
Case 3	$(-5,-3,\mathrm{\hspace{0.17em}0})$	$(5,\mathrm{\hspace{0.17em}3},\mathrm{\hspace{0.17em}0})$
Case 4	$(-5,-4,\mathrm{\hspace{0.17em}0})$	$(5,\mathrm{\hspace{0.17em}4},\mathrm{\hspace{0.17em}0})$
Case 5	$(-5,-5,\mathrm{\hspace{0.17em}0})$	$(5,\mathrm{\hspace{0.17em}5},\mathrm{\hspace{0.17em}0})$

Equations

$$f(\mathcal{E})=\left\{\begin{array}[]{ll}\displaystyle\frac{\rho_{1}}{{(2\pi\sigma^{2})}^{3/2}}(e^{\mathcal{E}/\sigma^{2}}-1)&(\mathcal{E}>0),\\ 0&(\mathcal{E}\leq 0),\end{array}\right.$$

$$\mathcal{E}=\Psi (r)-\frac{1}{2}({v_x}^2+{v_y}^2+{v_z}^2),$$

$$\Phi (r)=-\Psi (r)-\frac{GM}{r_t},$$

$$r_0=\sqrt{\frac{9{\sigma }^2}{4\pi G{\rho }_0}},$$

$${\nabla }^2{\Phi }_{nlm}(\text{r})=4\pi G{\rho }_{nlm}(\text{r}),$$

$$\int {\rho }_{nlm}(\text{r}){\Phi }_{n^{\prime }{l}^{\prime }{m}^{\prime }}(\text{r})d\text{r}={\delta }_{n{n}^{\prime }}{\delta }_{l{l}^{\prime }}{\delta }_{m{m}^{\prime }},$$

$$\rho (\text{r})=\sum _{n,l,m}{A}_{nlm}(t){\rho }_{nlm}(\text{r})$$

$$\Phi (\text{r})=\sum _{n,l,m}{A}_{nlm}(t){\Phi }_{nlm}(\text{r}),$$

$$\int \rho (\text{r}){\Phi }_{nlm}(\text{r})d\text{r}$$

$$A_{nlm}(t)$$

$$\text{a}(\text{r})=-\sum _{n,l,m}{A}_{nlm}(t)\nabla {\Phi }_{nlm}(\text{r}),$$

$$\text{a}({\text{r}}_{\mathrm{CM},i})=-\sum _k\frac{G{m}_k({\text{r}}_k-{\text{r}}_{\mathrm{CM},i})}{{|{\text{r}}_k-{\text{r}}_{\mathrm{CM},i}|}^3},$$

$${\text{a}}_{i\to i}(\text{r})=\sum _{n,l,m}{A}_{nlm,i}(t)\nabla {\Phi }_{nlm}(\text{r}-{\text{r}}_{\mathrm{CM},i}),$$

$${\text{a}}_{j\to i}(\text{r})=\sum _{n,l,m}{A}_{nlm,j}(t)\nabla {\Phi }_{nlm}(\text{r}-{\text{r}}_{\mathrm{CM},j}),$$

$${\text{a}}_{j\to i}(\text{r})=\sum _{n,l,m}{A^{(\mathrm{CM})}}_{nlm,j}(t)\nabla {\Phi }_{nlm}(\text{r}-{\text{r}}_{\mathrm{CM}}),$$

$${\text{a}}_i(\text{r})={\text{a}}_{i\to i}(\text{r})+{\text{a}}_{j\to i}(\text{r})(i\ne j).$$

$${\rho }_{nlm}(\text{r})=\frac{K_{nl}}{4\pi }\frac{a^{l+5}{r}^l}{{(a^2+r^2)}^{l+5/2}}{C}_n^{(l+1)}(\xi )\sqrt{4\pi }{Y}_{lm}(\theta ,\varphi ),$$

$${\Phi }_{nlm}(\text{r})=-\frac{a^{l+1}{r}^l}{{(a^2+r^2)}^{l+1/2}}{C}_n^{(l+1)}(\xi )\sqrt{4\pi }{Y}_{lm}(\theta ,\varphi ),$$

$$\xi =\frac{r^2-a^2}{r^2+a^2},$$

$$Y_{lm}(\theta ,\varphi )=\sqrt{\frac{2l+1}{4\pi }\frac{(l-m)!}{(l+m)!}}{P}_{lm}(\cos \theta )\exp (im\varphi ).$$

$$K_{nl}=4n(n+2l+2)+(2l+1)(2l+3).$$

Full text

galaxies: evolution — galaxies: kinematics and dynamics — methods: numerical

1 Introduction

$N$-body simulation is an indispensable method for studying astronomical objects whose constituents interact gravitationally with one another. In order to reveal detailed structures of a target object which is modeled with $N$ particles, we need to invest in it as many particles as possible. However, extremely large $N$-body simulations are prohibitively time-consuming because the number of force calculation per time step increases explosively as in $\mathcal{O}(N^{2})$ with increasing $N$, if the force exerted on a given particle is calculated by summing up the mutual gravitational forces over all other particles. This disadvantage has been ameliorated by adopting, e.g., a tree algorithm (Barnes & Hut, 1986) which reduces the number of force calculation to $\mathcal{O}(N\log N)$. Regarding this algorithm, a recently developed numerical library called FDPS has facilitated the implementation of a tree code running with a sufficiently high speed on a massively distributed-memory parallel computer (Iwasawa et al., 2016; Namekata et al., 2018), which will promote simulations with an ever-larger number of particles.

When an $N$-body method is applied to collisionless systems such as galaxies and clusters of galaxies, we should pay attention to the difference between a system composed of $N$ particles and its corresponding genuine collisionless system. In a strict sense, an $N$-body system in which the forces among particles are calculated, in essence, using pair-wise interactions is not collisionless. This is because the term “collisionless” means that particles interact with the mean force field generated by themselves, and so, the forces do not depend on the relative coordinates of pair-wise particles. However, by terminating an $N$-body simulation within a period much shorter than the two-body relaxation time even in that way of pair-wise force calculation, we are allowed to regard the system as collisionless in a practical sense. On the other hand, a mean-field approach ensures the collisionless nature at least formally. In reality, if a mean field is represented by a particle-based method, Poisson fluctuations induced by the finite number of particles cause collisional effects equivalent to two-body relaxation (Hernquist & Barnes, 1990). Nevertheless, such an approach is still desirable for collisionless stellar systems from some points of view. For example, once the force field is given, the orbits of all particles can be calculated independently, and thereby one $N$-body problem is reduced to $N$ one-body problems, as has been pointed out by Hernquist & Ostriker (1992). This approach is thus well suited to parallel computation as shown by Hernquist et al. (1995), so that extremely large-$N$ simulations will be made feasible. In addition, if the mean field is known at each time step, we can trace the orbits of particles starting from arbitrary positions in phase space, which enables us to reproduce phase space itself at a given time, as has been demonstrated by Hozumi (1997).

We can realize a mean-field approach to simulating collisionless stellar systems using expansion techniques for solving Poisson’s equation. Above all, a self-consistent field (SCF) method termed by Hernquist & Ostriker (1992), the idea of which dates back to that of Clutton-Brock (1972, 1973), is considerably useful. The essence of the SCF method consists in expanding the density and potential of a system in a set of basis functions. In this method, the expansion of the full spatial dependence makes the cpu cost proportional to particle number $N$. In addition, ideal load-balancing is easily achieved on a massively parallel machine owing to the perfect scalability inherent in the collisionless nature.

One of the disadvantages intrinsic to SCF methods is the limitation in the range of applicability because the expansion center is needed at every time step. By virtue of this inconvenience, SCF methods have been applied only to stellar dynamical problems in which systems of interest are isolated and fixed in space with the center of mass being immovable such as violent relaxation of stellar systems (Hozumi & Hernquist, 1995; Hozumi et al., 1996), halo evolution (Weinberg & Katz, 2002), and disk evolution (Hozumi & Hernquist, 2005; Hozumi, 2012). However, moving systems are ubiquitous in the real Universe, as gravitationally interacting events such as tidal encounter, collision and merging are commonly observed. In particular, merging is the fundamental process of hierarchical structure formation in the universe dominated by cold dark matter. Even after galaxies have fully grown up through hierarchical merging, they encounter or collide, and frequently result in merging to each other. Accordingly, many simulations have been devoted to reproducing the characteristic appearance of those tails, bridges, and antennae which are indicative of galaxy interactions (e.g., Toomre & Toomre, 1972; Lynds & Toomre, 1976; Theys & Spiegel, 1977; Farouki & Shapiro, 1982; Wallin & Stuart, 1992; Howard et al., 1993; Barnes, 1998; Naab & Burkert, 2003). Therefore, if SCF methods can be applied to merging systems, we will be able to understand in detail the features produced by merging processes and the properties of merged systems using the huge number of particles. As an instance, if the cusp found in the central region of a dark matter halo (Navarro et al., 1996, 1997) arises from the merging of clumps as suggested by Fukushige & Makino (1997), an SCF technique might help elucidate the origin of the power-law nature of the cusp that continues almost down to the center by assigning a tremendously large number of particles to each clump.

In conventional $N$-body methods, a softening length is introduced explicitly in order to avoid numerical divergences of the mutual force when two particles pass close to each other. On the other hand, an SCF method includes implicit softening for the force calculation in the sense that the force resolution is limited to a certain degree of precision because the expansion terms of the density and potential of a system are inevitably truncated at the finite numbers in a set of basis functions. However, the non-existence of an explicit scale length in interaction forces suggests that an SCF method may be able to describe the evolution of a stellar system under the pure Newtonian force law more faithfully than an explicitly softening-dependent method. In general, gravitational softening affects more severely the dynamics of so-called cold, rotation-dominated systems like galaxy disks (Miller, 1971, 1974; Earn & Sellwood, 1995) than that of hot systems like elliptical galaxies which are supported by velocity dispersion. In a sense, the orbital motion of two merging spherical galaxies is dynamically cold, even though each galaxy is a hot elliptical-like system. Consequently, the softened gravity might influence a merging process. In fact, we will uncover the effects of a softening length on the orbital phase of two merging systems.

In this paper, we present merging simulations of two spherical collisionless stellar systems on the basis of a mean-field approach which is realized using an SCF code developed here. In Section 2, models and initial settings employed are described. In Section 3, we explain how an SCF method is applied to merging simulations, along with the details of tree code simulations which are used not only for comparison but for making clear the effects of a softening length. Results are shown in Section 4. We discuss the effects of a softening length, computation time, and application problems of the SCF code for merging simulations. Conclusions are given in Section 5.

2 Models