Analytical model for the dynamical motion of the bulges of two interacting galaxies

TL;DR

This paper develops an analytical model to describe the dynamical motion of galaxy bulges during interactions, incorporating mass exchange and loss, and validates it against N-body simulations.

Contribution

It introduces a novel analytical approach to model galaxy interactions with mass exchange, validated by comparison with GADGET-2 simulation results.

Findings

Qualitative agreement with N-body simulations on galaxy separation evolution

Effective modeling of mass exchange and loss in galaxy interactions

Analytical model captures key dynamical features of interacting galaxies

Abstract

Two mathematical models of three bodies of variable masses are used to obtain a qualitative description of two interacting galaxies with mass exchange and mass loss. The reference system is centred on the largest body (the most massive galaxy), and the other two bodies are allowed to move around this one under the laws of gravity. The third body, which simulated the mass lost by the second galaxy in the form of a tail, increases its mass due to the mass lost by the second body and follows its trajectory. We are interested in knowing the time evolution of the separation of the two bulges of the interacting galaxies, and the parameters for the analytical models are obtained by running simulations with the GADGET-2 N-body code. The resulting behavior of this distance in our mathematical models is qualitatively in good agreement with that obtained by this code.

Figures

Equations

$$\frac{d(m_1{\mathbf{v}}_1)}{dt}=-G\frac{m_1{m}_2}{R_{12}^2}{\hat{R}}_{12}-G\frac{m_1{m}_3}{R_{13}^2}{\hat{R}}_{13},$$

$$\frac{d(m_2{\mathbf{v}}_2)}{dt}=+G\frac{m_2{m}_1}{R_{12}^2}{\hat{R}}_{12}-G\frac{m_2{m}_3}{R_{23}^2}{\hat{R}}_{23},$$

$$\frac{d(m_3{\mathbf{v}}_3)}{dt}=+G\frac{m_3{m}_1}{R_{13}^2}{\hat{R}}_{13}+G\frac{m_3{m}_2}{R_{23}^2}{\hat{R}}_{23},$$

$$\frac{d(m_2\mathbf{v})}{dt}=-G\frac{m_2{m}_1}{|\mathbf{x}{|}^{3/2}}\mathbf{x}+G\frac{m_2{m}_3}{|{\mathbf{x}}^{\prime }-\mathbf{x}{|}^{3/2}}({\mathbf{x}}^{\prime }-\mathbf{x}),$$

$$\frac{d(m_3{\mathbf{v}}^{\prime })}{dt}=-G\frac{m_3{m}_2}{|{\mathbf{x}}^{\prime }-\mathbf{x}{|}^{3/2}}({\mathbf{x}}^{\prime }-\mathbf{x})-G\frac{m_3{m}_1}{|{\mathbf{x}}^{\prime }{|}^{3/2}}{\mathbf{x}}^{\prime },$$

$$m_1(t)=m_{10}+(b_1+\beta )t,$$

$$m_2(t)=m_{20}+(b_2-\beta )t-m_3(t),$$

$$m_3(t)=m_{30}+m_{3f}(1-e^{-\gamma t}),$$

$$\frac{d(m_2\mathbf{v})}{dt}=-G\frac{m_2{m}_1}{|\mathbf{x}{|}^{3/2}}\mathbf{x}+G\frac{m_2{m}_3}{|{\mathbf{x}}^{\prime }-\mathbf{x}{|}^{3/2}}({\mathbf{x}}^{\prime }-\mathbf{x}),$$

$${\mathbf{x}}^{\prime }(t)=\mathbf{x}(t-t_r)+\alpha (t-t_r)\hat{\mathbf{x}}(t-t_r),$$

$$d(t)=|\mathbf{x}(t)|,$$

$$M_{200}=M_{bulge}+M_{halos}+M_{disk}.$$

$${\rho }_{halo,bulge}(r)=\frac{M_{halo,bulge}}{2\pi }\frac{a_{h,b}}{r(r+a_{h,b})},$$

$$J=\lambda \sqrt{2G{M}_{200}^3{r}_{200}/f_c},$$

$${\rho }_{\ast }(R,z)=\frac{M_{\ast }}{4\pi z_0{R}_d^2}sec{h}^2\left(\frac{z}{2z_0}\right)\exp \left(-\frac{R}{R_d}\right),$$

$${\Sigma }_{gas}=\frac{M_{gas}}{2\pi h_r^2}\exp \left(-\frac{r}{h_r}\right),$$

$$\frac{1}{{\rho }_{gas}}\frac{\partial P}{\partial z}+\frac{\partial \Phi }{\partial z}=0,$$

$${\mathbf{G}}_1:$$

$$c=13.0,v_{200}=220km/s,\lambda =0.033,M_{disk}/M_{200}=0.0,M_{bulge}/M_{200}=0.007,a_b/R_d=1.0,$$

$$N_{halo}=150,000,N_{disk}=0.0,N_{gas}=0.0,N_{bulge}=50,000,$$

$${\mathbf{G}}_2:$$

$$c=5.0,v_{200}=100km/s,\lambda =0.033,M_{disk}/M_{200}=0.041,M_{bulge}/M_{200}=0.01367,M_{gas}/M_{disk}=0.1$$

$$a_b/R_d=z_0/R_d=0.2,J_{disk}/J=0.041,$$

$$N_{halo}=300,000,N_{disk}=200,000,N_{gas}=20,000,N_{bulge}=10,000,$$

$$m_{10}=185.513000\cdot 10^{10}{\mathrm{h}}^{-1}{\mathrm{M}}_{\odot },m_{20}=17.4265690\cdot 10^{10}{\mathrm{h}}^{-1}{\mathrm{M}}_{\odot },m_{30}=m_{3f}=1.74265695\cdot 10^{10}{\mathrm{h}}^{-1}{\mathrm{M}}_{\odot },$$

$$\beta =1.30699262{\mathrm{f}}^{-1}{\mathrm{M}}_{\odot }{\mathrm{year}}^{-1},b_1=74.2051983{\mathrm{f}}^{-1}{\mathrm{M}}_{\odot }{\mathrm{year}}^{-1},b_2=21.7832112{\mathrm{f}}^{-1}{\mathrm{M}}_{\odot }{\mathrm{year}}^{-1},$$

$$\gamma =1.0{\mathrm{f}}^{-1}{\mathrm{h}\mathrm{Gyear}}^{-1},t_r=0.0{\mathrm{f}}^{-1}\mathrm{tGyear},|{\mathbf{x}}_0-{{\mathbf{x}}^{\prime }}_0|=24.3476219{\mathrm{h}}^{-1}\mathrm{kpc},$$

$$x_0=121.73811340332031{\mathrm{h}}^{-1}\mathrm{kpc},y_0=0.0{\mathrm{h}}^{-1}\mathrm{kpc},v_x(0)=0.0km/s,$$

$$v_y(0)=181.02100930177266km/s,v_x^{\prime }(0)=166.4443359375km/s,v_y^{\prime }(0)=221.70454793214924km/s,$$

$$m_{10}=185.513000\cdot 10^{10}{\mathrm{h}}^{-1}{\mathrm{M}}_{\odot },m_{20}=17.4265690\cdot 10^{10}{\mathrm{h}}^{-1}{\mathrm{M}}_{\odot },m_{30}=m_{3f}=1.74265695\cdot 10^{10}{\mathrm{h}}^{-1}{\mathrm{M}}_{\odot },$$

$$\beta =0.430285633{\mathrm{f}}^{-1}{\mathrm{M}}_{\odot }{\mathrm{year}}^{-1},b_1=137.417049{\mathrm{f}}^{-1}{\mathrm{M}}_{\odot }{\mathrm{year}}^{-1},b_2=21.5142822{\mathrm{f}}^{-1}{\mathrm{M}}_{\odot }{\mathrm{year}}^{-1},$$

$$\gamma =0.987654328{\mathrm{f}}^{-1}{\mathrm{h}\mathrm{Gyear}}^{-1},t_r=0.05{\mathrm{h}}^{-1}\mathrm{tGyear},|{\mathbf{x}}_0-{{\mathbf{x}}^{\prime }}_0|=24.3476219{\mathrm{h}}^{-1}\mathrm{kpc},$$

$$x_0=121.73811340332031{\mathrm{h}}^{-1}\mathrm{kpc},y_0=0.0{\mathrm{h}}^{-1}\mathrm{kpc},v_x(0)=0.0km/s,$$

$$v_y(0)=181.02100930177266km/s,v_x^{\prime }(0)=0.0km/s,v_y^{\prime }(0)=388.14888386964924km/s.$$

$$m_{10}=185.513000\cdot 10^{10}{\mathrm{h}}^{-1}{\mathrm{M}}_{\odot },m_{20}=17.4265690\cdot 10^{10}{\mathrm{h}}^{-1}{\mathrm{M}}_{\odot },m_{30}=0,m_{3f}=1.74265695\cdot 10^{10}{\mathrm{h}}^{-1}{\mathrm{M}}_{\odot },$$

$$\beta =0.430285633{\mathrm{f}}^{-1}{\mathrm{M}}_{\odot }{\mathrm{year}}^{-1},b_1=4.58056778{\mathrm{f}}^{-1}{\mathrm{M}}_{\odot }{\mathrm{year}}^{-1},b_2=21.5142822{\mathrm{f}}^{-1}{\mathrm{M}}_{\odot }{\mathrm{year}}^{-1},$$

$$\gamma =0.987654328{\mathrm{f}}^{-1}{\mathrm{h}\mathrm{Gyear}}^{-1},t_r=0.05{\mathrm{h}}^{-1}\mathrm{tGyear},{\alpha }_0=24.3476219{\mathrm{h}}^{-1}\mathrm{kpc},$$

$${\omega }_{\alpha }=0.987654328{\mathrm{f}}^{-1}{\mathrm{h}\mathrm{Gyear}}^{-1},x_0=121.73811340332031{\mathrm{h}}^{-1}\mathrm{kpc},y_0=0.0{\mathrm{h}}^{-1}\mathrm{kpc},$$

$$v_x(0)=0.0km/s,v_y(0)=181.02100930177266km/s.$$

Full text

Keywords: galaxies: kinematics and dynamics, galaxies: interactions, galaxies: bulges, galaxies: binary

PACS: 98.35.Gi, 98.35.Jk, 9865.At, 98.65.Fz

1 Introduction

Galaxies are one of the most important structures in our Universe, and studies of their formation and subsequent evolution continue to fill the research journals [White, 1978, Dekel, 1986]. All types of galaxy consist of the combination of an elliptical-like component (bulge) [Bouwens, 1999, Milosavljević, 2001, Bertin, 2000, Launhardt, 2002, Dwek, 1995, Xu, 2007, Salviander, 2007] and a disk-like component [van der Kruit, 2010, Bertin, 2000, Giallongo, 2000, Kennicutt, 1983, Kennicutt, 1987, Burkert, 1992], in different proportions (which can be zero), containing stars, gas, dust [Bertin, 2000, Kacprzak, 2011, Ivison, 2010, Daylan, 2016], plus a dark matter halo which usually contains the bulk of the mass. Visually the bulge component [Holtzman, 1998, Magorrian, 1998] defines the trajectory of the galaxy, usually coinciding with the centre of the dark matter halo.

Galaxies exist in binary systems and in small groups (less than 100 members) or large clusters of galaxies with increasingly complicated dynamics [Turner, 1976, Paul, 2017, Voit, 2005, Gallaghere, 1972]. When two galaxies interact, there is an exchange of mass [Gutowski, 1976], and usually one (or both) of the galaxies lose mass in the form of one [Kemp, 2016], or more usually, two tidal tails [Wilson, 1977, Vallée, 1976, Karachentsev, 2008]. Material in tidal tails may escape the system or fall back into the galaxy eventually [Hibbard, 1995, Kemp, 2016], or may be incorporated in tidal dwarf galaxies, galaxies newly formed from material in tidal tails [Duc, 2002, Kemp, 2016].The distribution of material around interacting galaxies and the forms of the galaxies themselves have been successfully modelled using numerical N-body codes since the works of [Toomre, 1972].

Analytically this phenomenon corresponds to a typical mass variation problem, which is called the Gyldén-Mestschersky problem in astronomy [Gyldén, 1884, Mestschersky, 1893, Mestschersky, 1902, Lovett, 1902, Bekov, 1988, Bekov, 1991],  which has had a long standing challenge for a complete analytical description [Berković, 1981, Bekov, 1989, Spivak, 2010, Sommerfeld, 1964, López, 2014].

There have been some works [Soares, 1996] with a semi-analytical description of the dynamics of two point galaxies san1,san2,bleh1,bitre,ezo but as far as we know, there have been few attempts to consider the mass variation problem in binary galaxies [Jeans, 1924, Jeans, 1928]. The analytical treatment of two interacting galaxies as a two-body problem cannot be considered satisfactorily since one needs to consider the dynamics of binary galaxies with mass loss or mass exchange between them.

Based on the observations made on previous studies [Spivak, 2010] and [López, 2013], we consider that Newton’s equations of motion for mass variation problems are still valid. In addition, we consider that the trajectory of a galaxy is defined mainly by its bulge (centre of the dark halo), and this is true even when two galaxies are interacting. We also consider that these bulges are treated as point bodies with variable mass, having one galaxy more massive than the other one. The less massive body loses mass (galaxy tails) which is going to be absorbed by a third body which also moves around the binary system, and it only increases its mass due to the mass loss from the second body. The parameters (mass variations) selected for our three bodies in the mass variation models will be taken from the N-body simulation program GADGET-2 [Spivak et al, 2005, Spivak, 2005, Springel, 1999] since the real evolution times of such systems are in Gyr. The principal point is to see if this very simple model of three bodies (galaxies) could qualitatively say something about the behavior of the very complicated real motion of two interacting galaxies, using parameters from the N-body simulations. The third body simulates a tidal tail, and if a second tail is present, it is effectively included in the mass of the first (most massive) body.

2 Analytical Models

Consider three bodies of initial masses $m_{1}(0)=m_{10}$, $m_{2}(0)=m_{20}$, and $m_{3}(0)=m_{30}$ under their mutual gravitational attraction. Then, the equations of motion of these bodies from an arbitrary inertial frame of reference are

[TABLE]

where $G$ is the gravitational constant ($G=6.674\times 10^{-11}$ m3kg*-1s-2*, [Rose, 1969]) ${\bf v}_{i}=d{\bf x}_{i}/dt$ for $i=1,2,3$ represent the velocities of the bodies, being ${\bf x}_{i}$ i=1,2,3 their positions, and $\widehat{R}_{ij}$ for $i\not=j$ are the unitary vectors $\widehat{R}_{ij}={\bf R}_{ij}/R_{ij}$, with $R_{ij}=|{\bf R}_{ij}|$, being ${\bf R}_{ij}={\bf x}_{i}-{\bf x}_{j}$.

2.1 Model-1

Moving our reference system to the body with mass $m_{1}$ (${\bf x}_{1}={\bf 0}$), the first equation is not considered, and defining ${\bf x}_{2}={\bf x}$ and ${\bf x}_{3}={\bf x}^{\prime}$, the following equations result,

[TABLE]

where the time depending mass will be taken as

[TABLE]

where $\beta$ represents the mass transfer between the bodies of masses $m_{1}$ and $m_{2}$, the parameters $b_{1}$ and $b_{2}$ represent the increases in the effective masses of the bodies $m_{1}$ and $m_{2}$ due to the increases of the local particle densities of the spheres in which we define the galaxy masses when using the GADGET-2 code, $m_{30}$ and $m_{3f}$ are the initial and final masses of the third body, and $\gamma$ is related to the rate of mass gain of the third body from the second one.

2.2 Model-2

In this model, we assume again that our reference system is located in the body with mass $m_{1}$, and that the trajectory of the body with mass $m_{3}$ follows the motion of the body with mass $m_{2}$ in a fixed specified way. Thus, equations (1a) and (1b) are not taken into account, and the motion of the system is reduced to the equation

[TABLE]

where the motion of the body with mass $m_{3}$ will be determined by the expression

[TABLE]

where $t_{r}$ is the retarded time which is the time delay for movement of the third body, after the second body has started to move (used to simulate the tail leaving the second galaxy). In this model, if $t_{r}=0$ the second and third bodies would be moving tangentially on the same radial line, separated by some defined distance, on their tangential motion on the same radial line. The function $\alpha(t-t_{r})$ is a function which takes into account the separation of the third body from the second one, $\alpha(t-t_{r})=\alpha_{0}(2-e^{-\omega_{a}(t-t_{r})})$.

Equations (2a), (2b), and (4) are solved numerically, with the functions $m_{i}(t)$ for $i=1,2,3$ determined qualitatively from the N-body GADGET-2 code, when two galaxies are interacting. The idea is to study the behavior of the separation between the body with mass $m_{1}$ and the body with mass $m_{2}$,

[TABLE]

and compare this parameter with the same parameter obtained using the GADGET-2 code.

3 Simulations of two interacting galaxies with GADGET-2

The GADGET-2 code is a program that models a galaxy and galaxy interactions. Each galaxy has a bulge, a disk, central black holes, gas, and dark matter halo. The description and capabilities of this program are given in [Spivak et al, 2005, Spivak, 2005, Springel, 1999]. If $r_{200}$ is the radius containing $75\%$ of the galaxy mass, where the density is 200 times the critical density and $v_{200}$ is the rotation velocity or velocity dispersion at this radius, then the total mass of the galaxy within this radius is $M_{200}=v_{200}^{2}r_{200}/G$, where $G$ is the gravitational constant.

[TABLE]

The halo and bulge densities, as a function of the radius ($r$) of the galaxy, follow the expression

[TABLE]

where the parameter $a_{h,b}$ is given in term of the concentration index $c$ as $a_{h}=r_{s}\sqrt{2[\ln(1+c)-c/(1+c)}$, and the halo scale length is defined as $r_{s}=r_{200}/c$ and the proportion bulge scale length per scale radius $f_{b}=a_{b}/R_{d}$, where $R_{d}$ is the scale length of the disk. The angular momentum of the halo is given by

[TABLE]

where $\lambda$ is the twist parameter, and the parameter $f_{c}$ is written in terms of the concentration index as $f_{c}=c[1-1/(1+c)^{2}-2\ln(1+c)/(1-c)]/(2[\ln(1+c)-c/(1+c)]^{2})$. The fractional angular momentum of the bulge is defined as $J_{b}=JM_{bulge}/M_{200}$. The disk (stars) density varies as

[TABLE]

where $R=\sqrt{x^{2}+y^{2}}$, $z_{0}$ is the parameter determining the thickness of the disk. The angular momentum of the disk is $J_{d}=j_{d}J$, with the free parameter $j_{d}$. Similarly, the gas in the disk is modelled as

[TABLE]

where $h_{r}$ is the scale length of the gas profile, and $M_{gas}+M_{*}=M_{disk}$. In addition, the vertical structure of gas in asymmetric galaxies is governed by the equation

[TABLE]

where $\Phi$ is the gravitational potential due to the total mass of the gas.

Our simulation of two interacting galaxies with the GADGET-2 code were performed with the following set of parameters for the galaxies $G_{1}$ and $G_{2}$ ($N_{i}$ being the number of elements of each part of the galaxy):

[TABLE]

where $R_{d}=3.7397$ h*-1*kpc.

[TABLE]

where $R_{d}=1.978$ h*-1*kpc.

Figure 1 shows the initial configuration of both galaxies with the regions (circles) representing spheres used to define the total masses of the galaxies in the planes x-y and x-z. The mass variation depends on the size of these circles. Figure 2 and 3 show the mass evolution of both galaxies during their interaction for several regions of different sizes (as the size increases, the number of elements increases and the mass increases, as shown initially in both figures) for the galaxies of mass $m_{1}$ and $m_{2}$, where the size of the regions are $50~{}$h*-1kpc and $50/2.2~{}$h-1kpc longer than the previous one, for G1 and G2 ($h$ is the normalization of the Hubble constant $H_{0}$, $h=H_{0}/100$ km s-1* Mpc*-1* and observationally, $h\sim 0.7$, and 1 kpc $\approx 206264806~{}$au $\approx 149.6\times 10^{6}~{}$km) The maxima on the mass of the second galaxy, shown in Figure 3, are due to excess or counting the elements when the two galaxies are colliding with each other, and the minima correspond when they are separated enough that this counting is well defined. Figure 4, green solid line, shows the evolution of the separation of the centre of the bulges of the two galaxies as a function of time, where we can see the shrinking damping oscillations of the system. To make the calculation of the distance between these two galaxies, we first located the centres of both galaxies (the points of highest density) and separate these centres a distance $d(0)=120~{}$h*-1*kpc in the x-direction. Then, we follow the evolution of the separation. This spatial scale covers the typical separations and spatial scales of filaments seen in interactiing galaxies like NGC 4038/4039 [Lahén, 2018] and merger remnants like NGC 7252 [Hibbard, 1995].

4 Analytical result and comparison

As we can see from Figures 1,2 and 3, we can consider the mass of the three bodies and the values of the parameters for our analytical Model-1 (without retarded time) as:

[TABLE]

where M⊙ represents the Solar Mass (M${}_{\odot}=1.989\times 10^{30}$kg) and the dimensionless parameter $\mathrm{f}=0.98$. For this same model but with retarded time ($t_{r}$), the parameters are:

[TABLE]

The parameters for our analytical Model-2 can be chosen as:

[TABLE]

Solving numerally the equations (2a), (2b), and (4), we obtain Figure 4, where the parameter (6) is plotted for the two mathematical models, Model-1 corresponds to the dotted green curve (no retarded time for the third body) and the dotted-dashed yellow curve (with retarded time $t_{r}=5~{}$h*-1Gyr for the third body), Model-2 corresponds to the dashed blue curve (with a retarded time for the second body of $t_{r}=5~{}$h-1*Gyr), and the result obtained with the GADGET-2 code (continuous red curve). As one can see, there is a qualitatively good agreement of the GADGET-2 code with both models, and Model-2 with Model-1 (with and without retarded time) have almost the same behavior, even though Model-2 is an oversimplifying approximation of Model-1.

Figure 5 and 6 show the motion of the second and third bodies (relative to the first body) on the plane of motion, and one can see that the third body tends to scape from the system (simulating the tail of the second Galaxy moving out of the binary system). Of course, one can not expect similar behavior since in Model-2 we have forced the third body to follow the second body at some retarded time ($t_{r}=5~{}h^{-1}$Gyr ). The small square mark on Figures 5 and 6 indicates the position of the second body when the third body is starting to move, due to retarded time.

5 Results and Comments

We have used two mathematical models for the restricted three-body problem of variable masses to obtain a qualitative description of two interacting galaxies with mass exchange and mass loss. The main point is the use of the third body as a ghost galaxy which simulated the mass lost by the second galaxy in the form of a tail, and this ghost galaxy increases its mass due to the mass lost by the second body. We were interested in predicting the time evolution of the distance separation of the two bulges of the interacting galaxies, using two separate models, one as two bodies with mass variation, the other with one body with mass variation, using parameters deduced by the use of the GADGET-2 N-body simulation code for two interacting galaxies. The general result we have obtained is that the time dependent distance between the bulges of the two galaxies, calculated with our mathematical models and deduced from the GADGET-code, agree well qualitatively using both our mathematical models. This agreement is better with our Model-1 (three-body problem, without retarded time) than with Model-2, but both have almost the same behavior. This implies that for a few orbits of two interacting galaxies (before they merge as a single bigger galaxy), we can use this very simple model to know qualitatively the trajectory of the bulges of the two galaxies. We need to point out that one needs just a single run on N-bodies code to obtain the necessary parameters for the analytical models and to know the dynamics of the budges in much more shorter time than the simulation code. As a final comment, we expect that including also the position depending of the masses of the bodies in our models, the results can be improve appreciably to match with simulation code.

Acknowledgements

We want to thank the INAOE institute in Tonantzintla, Puebla, México, to allow us to have access to their installations and for its financial support. In particular, we want to thank Dr. Ivânio Puerari and M.C. Diego Valencia Enríquez for their help teaching to one of us how to use GADGET-2 code.