Effects of galaxy--satellite interactions on bar formation

TL;DR

This study uses high-resolution simulations to show that low-mass satellite galaxy encounters can significantly delay or advance bar formation in Milky Way-like galaxies, with effects depending on satellite mass and orbit.

Contribution

It demonstrates how satellite interactions influence bar formation timing without affecting the final bar structure, highlighting the indirect role of low-mass satellites.

Findings

Satellite encounters can delay or advance bar formation by up to 1 Gyr.

Final bar properties are unaffected by satellite interactions.

Low-mass satellites below 10^9 M_sun have minimal impact unless very close.

Abstract

Aims. We aim to show how encounters with low-mass satellite galaxies may alter the bar formation in a Milky Way-like disc galaxy. Methods. We use high-resolution N-body simulations of a disc galaxy prone to mild bar instability. For realistic initial conditions of satellites, we take advantage of cosmological simulations of Milky Way-like dark matter haloes. Results. The satellites may have a significant impact on the time of bar formation. Some runs with satellites demonstrate a delay, while others show an advancement in bar formation compared to the isolated run, with such time differences reaching $\sim$ 1 Gyr. Meanwhile, the final bar configuration, including its very appearance and the bar characteristics such as the pattern speed and the exponential growth rate of its amplitude are independent of the number of encounters and their orbits. The contribution of satellites with masses below $10^9 M_{\odot}$ is insignificant, unless their pericentre distances are small. We suggest that the encounters act indirectly via inducing perturbations across the disc that evolve to delayed waves in the central part and interfere with an emerging seed bar. The predicted effect for the present-day host galaxy is expected to be even more significant at redshifts $z \gtrsim 0.5$.

Tables

Table 1. Table 1: Number statistics of Aq-D2 satellites.

Quantity	Value
$N_{\text{sub}}$ ($M_{\text{tid}}\ge {10}^6{M}_{\odot }$)	72,380
$f_{\text{sub}}$ ( $r_{\text{50}}$)	24.01 %
$f_{\text{cross}}$ ( $r_{\text{50}}$)	8.85 %
$N_{\text{cross}}$ ($>$ 108 $M_{\odot }$)	23
$N_{\text{cross}}$ ($>$ 3 $\times$ 108 $M_{\odot }$)	8
$N_{\text{cross}}$ ($>$ 5 $\times$ 108 $M_{\odot }$)	5
$N_{\text{cross}}$ ($>$ 109 $M_{\odot }$)	4

Table 2. Table 2: The characteristics of disk+bulge+halo for our host galaxy model.

Quantity	Value
${\Sigma }_{\text{sol}}$	50 $M_{\odot }$pc-2
$R_{\text{d}}$	2.9 kpc
z${}_{\text{d}}$	300 pc
$M_{\text{d}}$	4.2 $\times$ 1010 $M_{\odot }$
${\sigma }_R$	35 km s-1
${\sigma }_{R0}$	140 km s-1
$N_{\text{d}}$	6 million
$R_{\text{e}}$	0.64 kpc
${\sigma }_{\text{b}}$	272 km s-1
$M_{\text{b}}$	1.02 $\times$ 1010 $M_{\odot }$
$N_{\text{b}}$	1.5 million
$M_{200}$	1.29 $\times$ 1012 $M_{\odot }$
$r_{\text{s}}$	17.25 kpc
$r_{\text{200}}$	229.3 kpc
$N_{\text{h}}$	9.25 million

Table 3. Table 3: Summary of the runs

Model	$n_{\text{S}}$	${\Delta }_{\text{S}}$	$t_1$	$t_2$	$t_3$	${\Delta }_B$	${\Omega }_{\text{p}}$	${\omega }_{\text{I}}$
B-0	–	–	0.75	1.4	2.2	–	49	2.9
B-1	1	–	0.75	1.2	1.8	-200	48	3.0
B-1m	1	–	0.75	2	2.4	300	49	2.9
B-4	4	–	0.75	1.3	2.0	-300	49	2.9
B-4m	4	–	0.75	1.7	2.5	470	48	3.0
B-23	23	–	0.75	1.1	1.7	-310	48	2.9
B-1-110	1	+110	0.75	1.6	2.5	600	48	2.9
B-1-500	1	– 500	0.40	1.0	2.2	-20	48	2.7

Table 4. Table 4: Parameters of the bar (bar radius, pattern speed, corotation radius, ellipticity) in B-0 run at different moments in time.

Time [Gyr]	1.0	1.5	2.0	2.5	3.0	4.0
$R_{\text{b}}$ [kpc]	1.1	1.1	2.4	4.1	5.0	5.7
${\Omega }_{\text{p}}$ [km/s/kpc]	48.8	48.0	45.3	37.8	33.3	28.6
$R_{\text{C}}$ [kpc]	4.4	4.5	4.8	5.9	6.7	7.8
$\epsilon$	0.13	0.42	0.69	0.65	0.63	0.62

Equations

$$M=M_{\text{orig}}/f,$$

$$v_{\text{x,y,z}}=\frac{v_{\text{orig},x,y,z}}{f^{1/3}}\text{and}$$

$$x=\frac{x_{\text{orig}}}{f^{1/3}}y=\frac{y_{\text{orig}}}{f^{1/3}}z=\frac{z_{\text{orig}}}{f^{1/3}}.$$

$${\rho }_{\text{NFW}}(r)=\frac{{\rho }_{\text{s}}}{(r/r_{\text{s}})(1+r/r_{\text{s}}{)}^2}$$

$${\delta }_c=\frac{{\rho }_{\text{s}}}{{\rho }_{\text{crit}}}=\frac{200}{3}\frac{c^3}{\ln (1+c)-c/(1+c)}.$$

$${\rho }_{\text{b}}(r)={\rho }_{\text{b}}{\left(\frac{r}{R_e}\right)}^{-p}{\text{e}}^{-b(r/R_e{)}^{1/n}},$$

$${\sigma }_{\text{b}}\equiv {\left\{4\pi Gn{b}^{n(p-3)}\Gamma [n(3-p)]R_e^2{\rho }_{\text{b}}\right\}}^{1/2}.$$

$${\Omega }_{\text{p}}={\Omega }_{\text{pr}}(R)\equiv \Omega (R)-\frac{1}{2}\kappa (R).$$

$$Q=\frac{\kappa {\sigma }_R}{3.36G{\Sigma }_{\text{d}}},$$

$${\left(\frac{\Delta v}{v_{\text{c}}}\right)}_{R_{\text{d}}}\sim \frac{2G{M}_{\text{tid}}}{b^{2}V}\frac{R_{\text{d}}}{v_{\text{c}}}\sim \frac{G{M}_{\text{tid}}{R}_{\text{d}}}{b^2{v}_{\text{c}}^2}.$$

$$\mathbf{A}\mathbf{x}=\omega \mathbf{x},$$

$$A_m(t)=\sum _j{\mu }_j{\text{e}}^{-im{\theta }_j},$$

$$V_{k,m}(t)=\frac{1}{A_0}\sum _j{\mu }_j{\tilde{v}}_{k,j}{\text{e}}^{-im{\theta }_j},$$

$${\sigma }_k^2(t)=\frac{1}{A_0}\sum _j{\mu }_j{\tilde{v}}_{k,j}^2.$$

$$B(t)=1-I_{yy}/I_{xx},$$

$$\mathrm{d}\mathbf{v}(R,\theta )=-G\frac{M_{\text{s}}(t)(\mathbf{R}-{\mathbf{r}}_{\text{s}})}{|{\mathbf{r}}_{\text{s}}-\mathbf{R}{|}^3}\mathrm{d}t,\mathrm{d}\theta =\Omega (R)\mathrm{d}t,$$

$$v_{Ra}(R)=\frac{i}{\Delta}\Big{[}(\omega-m\Omega)\frac{{\rm d}\Phi_{a}}{{\rm d}R}-\frac{2m\Omega\Phi_{a}}{R}\Big{]}\ ,$$

$$v_{\theta a}(R)=-\frac{1}{\Delta}\Big{[}2B\frac{{\rm d}\Phi_{a}}{{\rm d}R}+\frac{m(\omega-m\Omega)\Phi_{a}}{R}\Big{]}\ ,$$

$$-i(\omega -m\Omega ){\Sigma }_{da}+\frac{1}{R}\frac{d}{dR}(R{v}_{Ra}{\Sigma }_0)+\frac{im{\Sigma }_0}{R}{v}_{\theta a}=0,$$

$$\Delta \equiv {\kappa }^2-(\omega -m\Omega {)}^2$$

$${\Phi }_a(R)=|{\Phi }_a|{\text{e}}^{iF(R)}=|{\Phi }_a|{\text{e}}^{i{\int }^{R}kdR},$$

$$k{v}_{Ra}{\Sigma }_0=(\omega -m\Omega ){\Sigma }_{da},$$

$$\begin{array}[]{rcl}F_{R}-F_{\Sigma}=\pi&:&R<R_{C}\ ,\\ F_{R}-F_{\Sigma}=0&:&R>R_{C}\ .\end{array}$$

$$v_{\theta a}{\Sigma }_0=\frac{-2iB}{\Delta }k{\Phi }_a,$$

$$F_{\theta }-F_{\Sigma }=-\pi /2,$$

$$v_{Ra}(R)=-\frac{2im}{\Delta }\frac{\Omega }{R}{\Phi }_a,$$

$$F_R-F_{\Sigma }=\pi /2,$$

$$v_{\theta a}(R)=-\frac{m}{\Delta }\frac{(\omega -m\Omega )}{R}{\Phi }_a,$$

$$\begin{array}[]{rcl}F_{\theta}-F_{\Sigma}=\pi&:&R<R_{C}\ ,\\ F_{\theta}-F_{\Sigma}=0&:&R>R_{C}\ .\end{array}$$

$$A_m(R,t)={\sum _j}^{\prime }{\mu }_j{\text{e}}^{-im{\theta }_j},$$

Full text

*Methods. *We use high-resolution N-body simulations of a disc galaxy prone to mild bar instability. For realistic initial conditions of satellites, we take advantage of cosmological simulations of Milky Way-like dark matter haloes.

*Results. *The satellites may have a significant impact on the time of bar formation. Some runs with satellites demonstrate a delay, while others show an advancement in bar formation compared to the isolated run, with such time differences reaching $\sim 1$ Gyr. Meanwhile, the final bar configuration, including its very appearance and the bar characteristics such as the pattern speed and the exponential growth rate of its amplitude are independent of the number of encounters and their orbits. The contribution of satellites with masses below $10^{9}M_{\odot}$ is insignificant, unless their pericentre distances are small. We suggest that the encounters act indirectly via inducing perturbations across the disc that evolve to delayed waves in the central part and interfere with an emerging seed bar. The predicted effect for the present-day host galaxy is expected to be even more significant at redshifts $z\gtrsim 0.5$.

1 Introduction

Structure formation in the Universe governed by the $\Lambda$ cold dark matter ($\Lambda$CDM) paradigm pursues a hierarchical model. Smaller DM (sub)haloes merge in order to build up larger objects (e.g. White & Rees 1978; Moore et al. 1999). Such mergers occur on different scales, ranging from minor (1:3 - 1:50 of total mass) to major (mass ratios $\gtrsim$ 1:3 of total galaxy mass) mergers. These events are believed to impact the dynamics of galaxies, and even change their morphology (e.g. Baugh, Cole & Frenk 1996; Borne et al. 1999; Hernández-Toledo et al. 2005; Lotz et al. 2010; Kannan et al. 2015). Current hydrodynamical cosmological simulations, which account for different small- and large-scale physical processes such as radiative cooling, feedbacks, star formation, magnetic fields and so on, are able to reproduce a spectrum of realistic galaxies with different morphologies, including ellipticals, barred, and unbarred spirals (Vogelsberger et al., 2014).

Observations of galaxies in the local and the high redshift Universe, have shown a high fraction of barred spirals (e.g. Van den Bergh et al. 1996; Abraham et al. 1999). The properties of the host disc galaxy is affected by the presence of a bar, which possesses a strong non-axisymmetric gravitational potential. For instance, the redistribution of stars/gas within the bar’s corotation radius, together with the evolution of the bulge, are strongly correlated with the bar (e.g. Hohl 1971; Kormendy 1982; Athanassoula 2005). Early N-body simulations were successful in generating spiral galaxies hosting bars (e.g. Miller, Prendergast & Quirk 1970; Ostriker & Peebles 1973; Athanassoula & Sellwood 1986).

According to observations, a higher fraction of barred spirals are found in denser environments, such as galaxy groups and clusters, than in isolation, suggesting the importance of tidal interactions for bar formation (Elmegreen, Elmegreen & Bellin, 1990). Noguchi (1987) pioneered numerical studies of bar formation in the presence of tidal interactions for initially stable galaxy models. In his runs, parabolic prograde planar very massive encounters with $M_{\textrm{sat}}/M_{\textrm{galaxy}}=1$ or 3 were able to only produce transient bars. Also, Walker et al. (1996) showed that a satellite with $M_{\textrm{sat}}/M_{\textrm{disc}}=0.1$ on a circular prograde orbit with a pericentre distance of 6$R_{\textrm{d}}$ and inclined by 30∘ with respect to the disc plane is capable of destroying the bar in the initially bar-unstable models (their galaxies transformed into Sa Hubble morphological type). In Moetazedian & Just (2016) we investigated the disc heating due to satellites infall in discs that were initially bar stable. In seven runs with different distributions of satellites extracted from zoom-in cosmological simulations of Milky Way-like hosts, no bar was induced, despite the mass ratios $0.003<M_{\textrm{sat}}/M_{\textrm{disc}}<4$. All these numerical experiments disfavour the tidally induced bar formation scenario in originally stable discs.

Dubinski et al. (2008) discussed that massive satellite galaxies crossing the inner region of stable host galaxies may excite a bar or a spiral structure in the disc through the process of Toomre’s swing amplification (Toomre, 1981, hereafter, T81), where the disc’s inner part becomes vulnerable to such amplifications.

Kazantzidis et al. (2008) analysed the impact of massive encounters on the disc galaxy evolution in the cosmological context. For $z=1$ epoch, the disc mass is expected to be lower, and the satellites distribution contains more massive satellites with smaller pericentre distances compared to the present day epoch $z=0$. The adopted disc mass and the exponential scale length are $M_{\textrm{d}}=3.53\times 10^{10}M_{\odot}$ and $R_{\textrm{d}}=2.82$ kpc, respectively. Also, the satellites were on eccentric orbits and possessed a mass range $0.21-0.57M_{\textrm{d}}$ with pericentre distances of $0.5-6.2R_{\textrm{d}}$. Unfortunately, the authors do not provide information concerning the stability of the isolated run. The models with satellites show growth of strong central bars, and the disc axisymmetry is not restored at late times. We note that the Toomre stability parameter $Q=2.2$ adopted for the isolated model is insufficient to avoid local non-axisymmetric (Polyachenko et al., 1997) and global instabilities such as bar instability as proved by our calculation of the isolated model with exactly the same Toomre parameter.

For the host galaxy in the present analysis we employ a stellar dynamical model with a cuspy bulge generated using the GalactICS code (Widrow et al., 2008), which was already studied in detail in our previous papers (Polyachenko et al., 2016a, b). It turns out that the process of bar formation is affected by randomly incoming waves towards the disc centre, delaying it for some time, and then facilitating the bar growth (Polyachenko, 2016).

The analysis with different disc mass favours a young bar hypothesis, according to which the bar instability is saturated only recently (Polyachenko et al., 2016b). Thus, the present paper is aimed at investigating the importance of satellite galaxy encounters on the time of bar formation in a Milky Way-like galaxy using high-resolution N-Body simulations, with disc parameters and satellite distribution closer to the present day epoch, in contrast to Kazantzidis et al. (2008). In order to have realistic initial conditions (ICs), we use the distribution of satellites with masses exceeding $10^{8}M_{\odot}$ extracted from cosmological simulations likely to host Milky Way-like galaxies (Springel et al., 2008; Moetazedian & Just, 2016).

In section 2 we discuss the cosmological ICs employed together with the characteristics of our desired host galaxy and the satellite galaxies. Section 3 gives a short overview of bar instability in the host galaxy. The spectral analysis employed here follows that described by Polyachenko (2005). The details of our simulations and the results obtained from our bar formation analysis are discussed in section 4. The relevance of the azimuthal phase between a weak bar and incoming satellites towards the bar formation epoch is investigated in section 5.1. Also, the results from the run with an early satellite encounter are presented in section 6. Section 7 contains comparisons with observations. We finalise the paper with a summary and discussion of the results in section 8.

2 N-body models

This section includes a description of the models recruited for setting up our desired Milky Way model, which consists of a disc, a bulge, and a DM halo together with the DM-only satellite galaxies.

2.1 Cosmological initial conditions

In order to account for a realistic image of satellite galaxies’ infall onto the Milky Way’s disc, the distribution of satellites was withdrawn from the Aquarius cosmological suite (Springel et al., 2008). This consists of a set of six realisations of DM haloes likely to host Milky Way-like galaxies and which have not experienced a recent major merger. For the purpose of this study the level 2 of the Aquarius-D simulation, hereafter Aq-D2, was employed. We have shown in Moetazedian & Just (2016) that Aq-D2 is a fair representative of a typical Milky Way-like system with contribution towards the vertical heating of the Galactic disc regarded as average compared to the rest of Aquarius simulations.

In our study we use the $z=0$ snapshot which corresponds to the present day distribution of DM substructures. The snapshot contains information such as position, velocity, maximum circular velocity $V_{\textrm{max}}$, the radial distance of this velocity $r_{\textrm{vmax}}$ and the tidal mass $M_{\textrm{tid}}$ for every substructure (satellite) within the simulation box. The parent host halo has an enclosed mass, $M_{\textrm{200}}=1.774\times 10^{12}M_{\odot}$, which corresponds to the mass within a sphere with radius $r_{\textrm{200}}=242.8$ kpc; this represents the radius at which the mean density of the DM halo is 200 times the critical density of the Universe ($\rho_{\textrm{crit}}$). Also the radius, where the halo encloses the mass with mean density 50 times $\rho_{\textrm{crit}}$ is denoted as $r_{\textrm{50}}$ and has a value of 425.7 kpc.

As mentioned earlier we are interested in satellites with a chance of passing close to the disc. The detailed procedure of identifying these candidates, is described in a previous paper (Moetazedian & Just, 2016). The statistics of the satellite distribution is shown in Tab. 1; $N_{\textrm{sub}}$ is the total number of satellites with $M_{\textrm{tid}}\geq 10^{6}M_{\odot}$ and $f_{\textrm{sub}}$ ($<r_{\textrm{50}}$) is the fraction of these objects within $r_{\textrm{50}}$. The percentage fraction of crossed satellites (i.e. objects that come closer than 25 kpc to the host halo’s centre during their 2 Gyr orbit) is represented by $f_{\textrm{cross}}$ ($<r_{\textrm{50}}$). The last four rows show the number of crossed satellites with masses $\geqslant$ 1, 3, 5 and $10\times 10^{8}M_{\odot}$. We argue in Moetazedian & Just (2016) that only satellites with $M_{\textrm{tid}}\geqslant 10^{8}M_{\odot}$ potentially contribute towards the heating of the Galactic disc. The Aq-D2 simulation has 23 satellites with $M_{\textrm{tid}}\geqslant 10^{8}M_{\odot}$ out of which four have masses larger than $10^{9}M_{\odot}$.

The host halo, in which these satellites would be inserted, has a slightly different enclosed mass ($M_{200}$). Therefore, we needed to rescale the phase-space properties of the satellites, together with their masses, following a recipe mentioned by Kannan et al. (2015), in order to have a physically sensible analysis. The scaling factor is defined as f=$M_{\textrm{200,Aq-D2}}/M_{\textrm{200}}$ (=1.37 in our case),

[TABLE]

[TABLE]

[TABLE]

The subscript “orig” represents the original non-scaled values. Early CDM N-body simulations have shown that the density profile of DM (sub)haloes could be well fitted via the known Navarro-Frank-White (NFW) profile (Navarro et al., 1997)

[TABLE]

and

[TABLE]

Here $\rho_{\textrm{s}}$ and $r_{\textrm{s}}$ are the scale density and the scale radius of the halo, respectively. The density contrast $\delta_{c}$ with respect to $\rho_{\textrm{crit}}$ is calculated with the concentration $c$ of the halo being the ratio of $r_{\textrm{200}}/r_{\textrm{s}}$.

As discussed in section 2.3, all the DM satellites also follow a NFW profile, which is tidally truncated at $r_{\textrm{tid}}$ according to the satellite mass $M_{\textrm{tid}}$.

2.2 The host galaxy

Our three-component model is adopted from Widrow et al. (2008); Kuijken & Dubinski (1995), and consists of a stellar disc, a bulge, and a DM halo. The complete list of characteristics for the disc, the bulge and the halo components are listed in Tab. 2. For a recent discussion of global parameters of the Milky Way see Bland-Hawthorn & Gerhard (2016).

The disc is exponential, with radial scale length $R_{\textrm{d}}=2.9$ kpc and truncation radius 15 kpc. The radial velocity dispersion $\sigma_{R}$ is approximately exponential with central value $\sigma_{R0}=140$ km s*-1* and radial scale length $R_{\sigma}=2R_{\textrm{d}}$. The solar neighbourhood location is at $R=8$ kpc (e.g. Gillessen et al. 2009; Reid et al. 2014).

The bulge component takes a density profile of the following form

[TABLE]

where $r$ is the spherical radius, corresponding to a Sérsic surface brightness profile with effective radius $R_{e}$. The scale density $\rho_{\textrm{b}}$, can be replaced by the bulge velocity scale

[TABLE]

Here, $\sigma^{2}_{\textrm{b}}$ corresponds to the depth of the gravitational potential associated with the bulge. In addition, we require $n=1.11788$ leading to $p\simeq 0.5$ parameter, in order to define the bulge density profile.

We have employed a truncated NFW profile for the halo, with the truncation radius at $r_{\textrm{200}}$. The bulge component dominates the inner rotation curve at $R\gtrsim 0.01$ kpc, despite the halo having a more cuspy profile.

The total circular velocity profile (solid curve) and the contribution from each component (dashed/dotted) are shown in the top panel (a) of Fig. 1. The bulge dominates at radii $R\lesssim 1.5$ kpc, while the halo takes over at $R>10$ kpc. At $R\approx 6$ kpc, where the contribution of the disc component reaches the maximum, the force from the halo is approximately two thirds of that from the disc in the Galactic plane.

Panel (b) presents the angular velocity profile $\Omega(R)$ in the equatorial plane. In addition, positions of the inner (ILR) and outer Lindblad resonances (OLR) are shown as $\Omega\pm\kappa/2$. For a given pattern speed $\Omega_{\textrm{p}}$, the ILR is calculated using

[TABLE]

The quantity, $\Omega_{\textrm{pr}}(R)$ is responsible for determining a precession rate of nearly circular orbits. Hence, it can be referred to as ‘precession’ curve, which diverges weakly as $R\to 0$ with $R^{-\alpha/2}$ and $\alpha\approx 0.5$.

The lower panel (c) of Fig. 1 demonstrates the Toomre $Q$ profile

[TABLE]

where $\Sigma_{\textrm{d}}$ represents the disc surface density. The solid line shows the profile which was initially set by the GalactICS code, while the dashed curve corresponds to the actual measured profile in our Galaxy run. The deviation for $2<R<12$ kpc is explained by larger actual radial dispersion. The initial minimum $Q_{\textrm{min}}=1.8$ resides at $R=5.9$ kpc. This minimum rises to $Q_{\textrm{min}}\approx 2.1$ in case of the actual run.

2.3 Satellites

In order to insert the satellites into our N-body simulations, each satellite needs to be generated as a distribution of particles following their corresponding NFW profile. The parameters $\rho_{\textrm{s}}$ and $r_{\textrm{s}}$ in Eq. 2.4 are determined by $r_{\textrm{vmax}}$ using $r_{\textrm{s}}=0.4623\,\,r_{\textrm{vmax}}$ and $V_{\textrm{max}}$ using the enclosed mass inside $r_{\textrm{vmax}}$. Next we derive the tidal radius $r_{\textrm{tid}}$ as cut-off radius from the cumulative mass profile reaching the satellite mass $M_{\textrm{tid}}$.

The distribution function introduced by Widrow et al. (2008) and implemented in Lora et al. (2013) was used for generating cuspy NFW profiles. Such a profile has an infinite cumulative mass at $r\rightarrow\infty$; therefore, we use the satellites’ $r_{\textrm{tid}}$ as the cut-off radius. In this work, all satellites are represented using 50,000 particles.

The satellites can be ranged according to their ability to excite density waves in the disc of the host galaxy. In order to quantify it, we shall introduce a ‘tidal impact’ parameter, which is a normalised tidal velocity perturbation. For a crude estimate, we use the formulae for tidal shocks (Binney & Tremaine, 2008, sect. 8.2.1), calculating the ratio of the perturbation to the circular velocity at some typical radius, for example, $v_{\textrm{c}}(R_{\textrm{d}})$:

[TABLE]

Here $b$ is the pericentre distance and we assumed that velocities at the time of encounters $V$ are similar and approximated as $V\simeq 2v_{\textrm{c}}$. Fig. 2 shows the reverse cumulative histogram of the number of encounters with the tidal impact parameter (2.10) greater than a given one, for all encounters with $M_{\textrm{tid}}>10^{8}M_{\odot}$ among seven runs discussed in Moetazedian & Just (2016). The Aq-D2 run can be regarded as a natural choice, since it contains at least a couple of satellites with tidal impact parameters at the high end. The marked red arrow corresponds to the Aq-D2 satellite with the highest parameter (referred to hereafter as ‘primary’), coming as close as $\sim 4$ kpc to the centre, and having the second highest mass in Aq-D2. The primary satellite is expected to have a larger tidal impact on the disc, by a factor of $\sim 2$, compared to the most massive satellite.

The host halo has a different enclosed mass and concentration compared to the Aq-D2 main halo. This means the orbits of Aq-D2 satellites would differ if inserted into this model. In order to gain fair and physically sensible results, we decided to scale $M_{\textrm{tid}}$, positions and velocities of Aq-D2 satellites using the factor $f=1.37,$ which corresponds to the ratio of Aq-D2 and the targeted host halo mass. Simulations with these initial conditions are called ‘fully-rescaled’. For comparison we also performed simulations where only the masses of the satellites are rescaled, but not their initial orbital positions and velocities, which we call ‘mass rescaled’ (see Tab. 3).

3 Bar instability in the host galaxy

The isolated galaxy is prone to instability leading to the formation of a bar. Bar properties on low amplitudes, such as a pattern speed and an amplitude (exponential) growth rate, can be obtained as global solutions from matrix equations (e.g. Kalnajs, 1971, 1977), describing a disc with stellar orbits of different eccentricities. It opposes tightly wound local WKB solutions of Lin–Shu–Kalnajs obtained in the epicyclic approximation.

The applicability of the razor-thin disc model for the description of real discs with finite thickness and a density cusp in the centre has been studied in detail by Polyachenko et al. (2016a). The main idea is to consider an effective rotation curve, which takes into account the $z$-dependence of the radial force on the stars that elevate above the equatorial plane. The effective rotation curve lacks cuspy features and possesses a maximum, which allows for a bar in a wide range of pattern speeds.

Here we use a matrix method by Polyachenko (2005) that has a form of the linear matrix equation,

[TABLE]

which allows us to find unstable modes effectively without prior information on the localisation of modes.

A serious flaw in all matrix methods is the inability to calculate unstable global modes in a live halo and bulge. Although the substitution of a live halo with its rigid counterpart gives nearly the same pattern speeds as in the fully live models, the growth rates are smaller by a factor of two or even more. The theoretical bar pattern speed estimate for our host galaxy is $\Omega_{\textrm{b}}\approx 48.34$ km s*-1* kpc*-1* and for the growth rate we obtained $\omega_{\textrm{I}}\approx 1.23$ Gyr*-1*, so the expected value for the fully live models is 2…3 Gyr*-1*.

4 N-body simulations

The initial conditions of stars for single mass simulations have been generated by the ‘GalactICS’ code provided by Widrow et al. (2008). In our previous N-body simulations (Polyachenko et al., 2016a) we used three different codes (SUPERBOX--10, bonsai2, ber-gal0) and demonstrated the practical independence of the results on the chosen code. The simulations in this work are restricted to the bonsai2 tree-code and ber-gal0333ftp://ftp.mao.kiev.ua/pub/users/berczik/ber-gal0/ (Zinchenko et al., 2015) routines as an auxiliary tool to analyse the snapshot data. The modified version of the recently developed N-body Tree–GPU code implementation bonsai2 (Bédorf et al., 2012a, b) includes expansion for force computation up to quadrupole order. The opening angle used had a value of $\theta=0.5$, accompanied by individual gravitational softening of 10 pc. bonsai2 employs the leap-frog integration scheme with a fixed time–step $\Delta t=0.2$ Myr.

The current set of simulations were carried out with the GPU version of the code using the ARI GPU cluster kepler and also the GPU cluster MilkyWay specially dedicated to the SFB 881 (“The Milky Way System”), located in the Jülich Supercomputing Centre in Germany.

Tab. 3 contains the summary of our runs. The capital letter ‘B’ denotes the numerical bonsai2 Tree code. It is worthwhile to mention, the Milky Way halo components have multi-mass halo particles in order to achieve a better resolution in the bar region. For this, we modified the GalactICS code utilising the same strategy as described in Dubinski et al. (2009). In the region between 0.1 and 1 kpc, the number density ratio of our multi-mass and single mass runs varies from 10 to 100, thus the effective numerical resolution there is enhanced by this factor.

The subscript ‘m’ corresponds to the simulations, which are only mass-rescaled using Eq. 2.1; while the remaining simulations are fully-rescaled. In the case of simulations marked using ‘B-1’, we use the primary satellite from Aq-D2 simulation with rescaled mass $M_{\textrm{tid}}=2.45\times 10^{9}M_{\odot}$. For runs marked as ‘B-4’, the four most massive Aq-D2 crossed satellites with $M_{\textrm{tid}}>10^{9}M_{\odot}$ were employed.

The orbital positions and velocities of the satellites are determined by their density centres. Fig. 3 shows the galactocentric orbital distance for Aq-D2’s primary satellite from our five runs. There exists a 110 Myr difference between the first pericentre passage of B-1 and B-1m runs. The pericentre distances have the values of 2.6 and 3.5 kpc, occurring at 0.85 and 0.96 Gyr, respectively. Also, the orbits from the two simulations with $+110$ and $-500$ Myr time shift are shown.

4.1 Isolated simulation B-0 of the host galaxy

In the absence of the satellites, the host galaxy appears unstable as manifested in the formation of a bar. In Fig. 4 one can see snapshots at times 2, 2.5, 3, and 4 Gyr of the bar oriented along the $x$-axis. At earlier times, the asymmetry of the density distribution is weak, because of a delay in bar formation (e.g. Polyachenko et al. 2016a; Polyachenko 2016). Moreover, there are no noticeable spirals throughout the simulation.

The bar is a wave, and its growth can be investigated using relative amplitudes of perturbed quantities, such as the surface density $A_{2}/A_{0}$,

[TABLE]

(here $\mu_{j}$ and $\theta_{j}$ are mass and polar angle of star $j$; $j$ spans particles within some fixed radius, for example, the radial scale length $R_{\textrm{d}}$), or velocity components $U_{k,m}$. The latter are defined as the perturbed velocity amplitudes,

[TABLE]

normalized to the average velocity dispersion $\sigma_{k}(t)$ calculated by

[TABLE]

The cylindrical components $R$, $\theta,$ and $z$ correspond to the quantity $k,$ and $\tilde{v}_{k,j}$ is the residual velocity component of the particle $j$.

Panel (a) in Fig. 5 shows the relative bisymmetric amplitudes $A_{2}/A_{0}$, and $U_{k,2}$ for the B-0 run, representing the isolated galaxy. All the curves, except $U_{z,2}$, show a typical phase of low amplitude fluctuations, followed by a jump. The evolution then enters a phase characterised by a typical exponential growth, and finally reaching a plateau. In Tab. 3, $t_{1}$ denotes the jump, $t_{2}$ corresponds to the start of the exponential growth and $t_{3}$ marks the end of the exponential growth. The simultaneous growth of the perturbed density and velocity components in the plane is expected for unstable density waves, in particular for bars. However, the vertical velocity component $U_{z,2}$ does not respond to the bar formation, and we exclude it from our analysis.

The calculation of components of the inertia ellipsoid can also be used for the analysis of the bar growth, and to obtain the pattern speed of the bar. Slopes of the bar amplitudes and the bar strength,

[TABLE]

give estimates for the growth rates, which appear to be very close to one another (and $\simeq 3$ Gyr*-1* as predicted from the rigid halo analysis).

The pattern speed is obtained from an angle of the rotation of the ellipsoid’s main axes (e.g. Wu et al. 2016). In panel (b) of Fig. 5 we plotted the pattern speed of the bar $\Omega_{\textrm{p}}$ determined by an average centred finite difference over two time-intervals of $\pm 10$ Myr before and after the evaluated time (dots) and the pattern speed after filtering out the highest frequencies. No definite structure for the pattern speed value can be found before the jump, $t<t_{1}$. Just after the jump, the centred finite difference shows much less scatter, and the filtered $\Omega_{\textrm{p}}$ varies only slowly until $t=1.9$ Gyr. Due to a bar slowdown, the pattern speed gradually decreases to 28.6 km s*-1kpc-1* at $t=4$ Gyr.

Fitting ellipses to the isophotes is one of the possibilities to quantify the bar shape (Martinez-Valpuesta et al., 2006). We define a bar radius, $R_{\textrm{b}}$, as a radius at which the ellipticity $\varepsilon(r)\equiv 1-b_{\textrm{e}}/a_{\textrm{e}}$ ($a_{\textrm{e}}$ and $b_{\textrm{e}}$ are the major and minor semi-axes of ellipses) declines by $\sim 15$% from its maximal value. The obtained radii and ellipticities are given in Tab. 4. For $t=1.0$ Gyr, the bar radius is approximately 1.1 kpc, and $\varepsilon\simeq 0.13$. At the end of the exponential growth, $t\approx 2.5$, the bar radius is 4.1 kpc, with an ellipticity of $\varepsilon\simeq 0.65$.

With the determined $\Omega_{\textrm{p}}\approx 49$ km s*-1* kpc*-1* at $T=1$ Gyr, the corotation resonance occurs at $R_{\textrm{C}}\approx 4.4$ kpc. The ILR and OLR are located at 0.7 and 7.9 kpc, respectively. As the bar grows between $t_{2}$ and $t_{3}$, the ratio ${\cal R}\equiv R_{\textrm{C}}/R_{\textrm{b}}$ attains the value of 1.4 and remains roughly constant until the end of the simulation despite the significant slow-down of the bar in the saturation phase.

As can be seen from Fig. 5a, the bar instability in the real N-body finite thickness disc does not manifest itself as an exactly exponential amplitude growth from the early beginning of the simulation. For a long time, the amplitude is on the noise level. The start of a mild growth is possible only after a sufficiently strong wave appears in the disc at some moment $t=t_{1}$, and provokes a delayed perturbation in the centre (see also top panels in Fig. 11). The emergence of such a wave perturbation is a random event, the probability of which is smaller for less unstable discs (Polyachenko, 2016).

The delayed wave is a bar-like perturbation which rotates with a well-defined pattern speed (Fig. 5b), suggesting that a seed bar is already formed at $t_{1}$. However, the amplitude growth now is irregular and slower than exponential, until the wave amplitude $A_{2}$ reaches 1 or 2 per cent of the axisymmetric background $A_{0}$ ($t=t_{2}$). The exponential growth with a growth rate $\omega_{\textrm{I}}\approx 3$ Gyr*-1* lasts until $t=t_{3}$, when the instability saturates.

4.2 Simulations with satellites

Possible effects from satellite interactions were studied through five runs, B-1, B-1m, B-4, B-4m and B-23. The numbers in the run codes reflect number of the satellites, while ‘m’ reflects the runs in which only the mass of the satellites was rescaled.

In Fig. 6 we compare the bar amplitudes $A_{2}/A_{0}$ of these models with the one calculated for B-0 (solid red). The jumps in all curves occur at the same time, $t_{1}=0.75$ Gyr, meaning that it is independent of the satellites, but rather intrinsic to the disc itself. The subsequent behaviour is different: mass rescaled runs B-1m and B-4m show a delay in bar formation, while in the fully-rescaled runs, B-1, B-4, and B-23, the bar forms earlier compared to B-0. The time difference in the bar formation $\Delta_{B}$ is given in Tab. 3.

In the mass-rescaled runs, the initial positions of the satellites are further away from the host galaxy than in the fully-rescaled runs, so the interaction with the disc peaks later in the mass rescaled runs. In particular, the minimum distance in the B-1m run is at 0.96 Gyr, while in the B-1 run it occurs at 0.85 Gyr. The delay $\Delta_{B}$ is larger in the B-4m run where four satellites are used.

The fully-rescaled models show contrary effect on the bar formation, compared to the mass-rescaled simulations. In the B-1 run the satellite advances the bar formation by 200 Myr (negative delay), and the higher the number of satellites, the larger the impact. In the B-4 run where the four largest satellites are used, the time difference is 300 Myr. A comparison with the B-23 run shows that the most massive satellites determine mainly the time difference – it is only 10 Myr larger than in the B-4 run.

In section 5 we investigate in more details the true character of the observed delay/advancement.

5 Phase relation between bar and satellite perturbations

A conjecture we wish to explore is focused on the measured mutual phase of an initial bar and perturbations induced by the massive satellites, and the importance towards advancement or delay in bar formation. Within the course of this section we inspect phases calculated using the surface density and velocity perturbations, assuming a simplified description of cold fluid discs.

5.1 Perturbations from a satellite

In this subsection we concentrate on the velocity perturbations in a cold fluid disc as a result of a satellite passage. In our runs, the typical encounter speed is $V\approx 400$ km s*-1*; for a galaxy of size $d\sim 20$ kpc the encounter lasts $d/V\approx 50$ Myr, which is a fraction of the rotation period of the bar, $\approx 130$ Myr, measured at the beginning of its formation.

To obtain velocity perturbations, one can integrate the acceleration of fluid particles moving on exactly circular orbits:

[TABLE]

where $M_{\textrm{s}}(t)$ and ${\mathbf{r}}_{\textrm{s}}(t)$ are known functions for the satellite mass and its distance from the disc centre, ${\mathbf{R}}$ is the radius vector of a fluid particle. We take only the in-plane acceleration into account and disregard vertical perturbations here.

The phase patterns, which are irregular at the beginning of the integration, become steady and regular as the satellite approaches the disc. In Fig. 7, we show amplitudes and phases of the bisymmetric $m=2$ radial and azimuthal velocity perturbations during the encounter of the primary satellite for the period between 0.82 and 0.87 Gyr. Loci of the maxima of $V_{R,2}$ and $V_{\theta,2}$ are practically unchanged from $t=0.7$ to 0.84 Gyr. After the fly-by, the maxima lines begin to wind up, and eventually turn into a tightly wound spiral. The amplitudes of velocity components grow up to $t=0.86$ Gyr, then remain practically unchanged.

5.2 Phase relations in the cold fluid disc

Equations of fluid dynamics are employed frequently to analyse phenomena taking place in stellar discs (e.g. Bertin, 2014). Indeed, when the stellar disc is cold enough, mean radial and azimuthal velocities can be obtained from linearised Euler equations for zero-pressure discs, because a cold stellar disc is dynamically equivalent to a fluid disc with zero pressure. If the perturbation has $m$-fold rotation symmetry and $\propto\textrm{e}^{i(m\theta-\omega t)}$, the perturbed surface density $\Sigma_{da}$ and average velocity components ($v_{\theta a},v_{Ra}$) obey the following relations (Binney & Tremaine, 2008):

[TABLE]

[TABLE]

[TABLE]

where $B\equiv-\kappa^{2}/4\Omega$ is the Oort constant; $\Sigma_{0}$ is the disc surface density and $\Phi_{a}$ is the potential perturbation;

[TABLE]

is positive between the Lindblad resonances.

In this analysis, we shall distinguish two extreme cases: tightly wound spirals and the bar. In the tightly wound approximation, a perturbed quantity (e.g. the potential) can be written in the form

[TABLE]

with phase function $F(R)$, where $k=dF(R)/dR$ and $|kR|\gg 1$. From eq. 5.4 one has

[TABLE]

that is, surface density of the wave is in phase with the radial velocity outside corotation, and in anti-phase inside corotation:

[TABLE]

The equation for the azimuthal velocity (eq. 5.3) simplifies to

[TABLE]

and taking into account that $\Phi_{a}=-2\pi G\Sigma_{da}/|k|$,

[TABLE]

meaning that the surface density lags behind the azimuthal velocity by $\pi/4$ for an $m=2$ perturbation.

In the bar region we assume that the radial derivative is negligible compared to the other term in square brackets in (5.2, 5.3). For the radial velocity,

[TABLE]

that is,

[TABLE]

meaning that radial velocity lags behind the surface density by $\pi/4$ ($m=2$). For the azimuthal velocity,

[TABLE]

that is,

[TABLE]

5.3 Phase relations in the stellar discs

To compare the phases of surface density and velocity perturbations in the stellar discs, we calculate a relative surface density, $A_{2}(R,t)/A_{0}(R,t)$,

[TABLE]

and velocity components, $V_{k,2}(R,t)/A_{0}(R,t)$,

[TABLE]

In both expressions, the summations are taken over disc particles within a ring $\Delta R$ near $R$.

Fig. 8 shows the loci of maxima $\theta(R)=-F/m$ for the B-0 run without satellites. Note that the maxima of $V_{\theta,2}$ are shifted by $\pi/2$ from the surface density maxima, and the maxima of $V_{R,2}$ lag by $\pi/4$ from the surface density maxima, in accordance with the phase relations for a bar in the cold fluid disc derived above.

The comparison with Fig. 7 shows that the wave velocity components are approximately in phase with disc perturbations caused by a satellite when the satellite approaches the disc centre, but become out of phase after $t=0.85$ Gyr.

Fig. 9 is similar to Fig. 8, but for the primary satellite in the mass scaled models. The orientation of the bar is roughly perpendicular to that of the bar in Fig. 8 at the corresponding moments of time (time difference is $\sim 110$ Myr). So in this case the bar is out of phase when the satellite approaches the pericentre.

The derived difference in phases during the satellite approach is essential, but seems to be only part of a more complex story. It follows from Fig. 6, where one can see a rather complicated bar amplitudes’ behaviour. In Fig. 10 we show a zoom-in of the early growth phase for the B-0, B-1, and B-1m runs. In particular, we stress the following features.

In the fully-rescaled B-1 run, the maximum occurs at $t=0.97$ Gyr, that is, 120 Myr (almost one period of bar rotation) after the pericentre passage; the 0.97 Gyr maximum is followed by minima at 1.06 Gyr and 1.16 Gyr. In the mass-rescaled model B-1m, a deep minimum occurs at $t=1.01$ Gyr, that is, 50 Myr after the pericentre passage; the minimum is followed by a large maximum exceeding the amplitude in the B-0 run at 1.11 Gyr followed by a nearly constant behaviour until 1.5 Gyr. These features cannot be understood from simple phase arguments as discussed above.

5.4 Interference of bar and satellite perturbations

The interaction between the central bar-mode and a perturbation from the satellite can be studied using density surface maps showing the total bisymmetric ($m=2$) amplitude in the $(R,t)$-plane. Left panels in Fig. 11 show $\log|A_{2}/A_{0}|$ for the isolated run B-0, as well as for one satellite runs B-1 (fully-rescaled) and B-1m (mass-rescaled). The lower left panel contains the map for not-yet-mentioned B-1-110 run, which differs from B-1 by a 110 Myr time delay of the satellite encounter. The lag is calculated to mimic a passage of the satellite near the pericentre of the mass-rescaled B-1m run (see Fig. 3).

The descending ‘tongues’ observed in the panels are nothing but interference patterns of two waves. Even in the upper figure for the isolated B-0 run, it is the interference between a small bar-like perturbation mostly residing in the central part, and random perturbations of the disc at $R\sim 1...2\,R_{\textrm{d}}$. These perturbations may act in accordance and accelerate bar formation, or, on the contrary, destroy one another.

The random peripheral perturbations are repeated in the lower panels. In addition, they contain imprints from the satellite for the different cases B-1, B-1m, B-1-110. In Fig. 14 one can see for case B-1 the onset of the disc perturbation after the encounter (0.87 and 0.89 Gyr), which clearly has a bar-like global shape and obeys phase relations obtained for bars in Section 5.2. This disc perturbation evolves so that its shape winds up gradually to the tightly wound spiral (last row of Fig. 14; the subsequent evolution is not shown here), but also it causes a disc response in the central region similar to those presented in T81; we refer to Figs. 1, 2 therein. Recall that Toomre analysed responses from weak corotating (Fig.1) and immobile (Fig.2) transient imposed dipole forces in a stable Mestel disc. The disc shows strong spiral responses developing from the centre outward, and sheared with a speed comparable with material arm shear, $\Omega(R)$. This is a manifestation of ‘a strong cooperative effect’ called swing amplification. For this mechanism to be efficient, the parameter

[TABLE]

should not exceed 3, which is the case in our model in the range $0.6<R<5.5$ kpc, with a minimum $X\simeq 2$ at $R\simeq 2.3$ kpc.

Much like in Toomre simulations, we expect similar responses to the disc external perturbations imposed by the satellite. Such responses grow on the time-scale of one period in the centre and interact with a small bar, strengthening or weakening the latter.

Right panels of Fig. 11 show $\textrm{Re}\,A_{2}/A_{0}$, that is, a joint $m=2$-amplitude along the $x$-axis, for the same runs B-0, B-1, B-1m, B-1-110 from top to bottom. In the upper panel, one can see red and blue spots standing for maxima and minima and showing emergence of the bar, starting from $t\sim 0.75$ Gyr ($R<2$ kpc). These spots extend vertically in time as the bar length grows. Also one can note smaller slightly shifted spots above the bars. These are interference patterns between the shearing external disc response to the encounter and the spiral arm adjacent to the rotating bar in the centre.

The interference spots are localised roughly along the straight lines with inclination growing in time. Using these lines one can infer periods of the interference pattern as a function of radius, which roughly obey the linear law, $T_{\textrm{i}}=aR+b$. Since frequency difference between the bar and the spiral response $\Omega_{\textrm{b}}-\Omega_{\textrm{s}}(R)=2\pi/T_{\textrm{i}}$, one can calculate the local pattern speed $\Omega_{\textrm{s}}(R)$ for wave patterns due to the encounter. This is shown by red circles in Fig. 12, along with $\Omega$ and $\Omega-\kappa/2$ curves.

Local pattern speeds $\Omega_{\textrm{s}}$ of the wave patterns can be explained using the simplest fluid WKB dispersion relation. Each wave pattern evolves according to the relation,

[TABLE]

meaning that $R=R(k,\omega)$ is a function of $k$; $\omega\equiv m\Omega_{\textrm{s}}$ is a parameter. For very open spiral555Application of WKB theory should be treated here as extrapolation only., $k$ is close to zero, so $\Omega_{\textrm{s}}=\Omega(R)-\kappa(R)/m$, and $R$ is close to the ILR, $R_{\textrm{ILR}}$. As the local wave pattern evolves, $(\omega-m\Omega)^{2}$ attains its minimum at some $k=k_{*}$, so $R$ attains its maximum (Toomre, 1969). Then $k_{*}$ can be estimated from the dispersion relation for fluid discs:

[TABLE]

where $k$ is the wavenumber and $c$ is a sound speed, which we use interchangeably with the radial velocity dispersion. Thus

[TABLE]

and substitution to the WKB relation gives

[TABLE]

In the last expression, $Q$ denotes Toomre stability parameter for fluid discs,

[TABLE]

As is seen from Fig. 12, $\Omega_{*}$ and $\Omega_{\textrm{s}}$ agree relatively well, which proves that the picture presented here is an interference pattern between the bar and the density wave. We note that due to large $Q$ ($>2$), radial excursions of the wave patterns are small, and are localised close to their ILR.

The interference of the bar and the delayed response to the disc perturbation in the B-1 run leads to a bar enhancement. This is seen from comparison of the right panels in the first and second rows of Fig. 11. In particular, a red bar spot immediately after $t=1$ Gyr in the B-1 run is definitely stronger than its counterpart in the B-0 run.

Contrary to simulations performed in T81, our external perturbation does not appear for only a short time. Thus it provides a number of delayed responses in the centre which can also either enhance or destroy the bar, unless the external perturbation becomes tightly wound. In case of B-1, it also works for enhancement.

In the last two rows, we see the opposite behaviour. In B-1m, the red bar spot immediately after $t=1$ Gyr is definitely weaker than in B-0 run, and it is seen that the delayed response is not in phase with the original bar. Moreover, it seems that the established perturbation also suffers from the subsequent delayed response waves.

We also note powerful random waves seen in B-0 panels between 1.05 and 1.3 Gyr, and 1.15 and 1.4 Gyr. They also came into play and decrease amplitudes in B-1m and B-1-110 runs.

6 The early encounter

In the previous section we argue that satellite encounters may cause a delay or advancement in bar formation after the small bar is already formed in the centre ($t>t_{1}$), but its amplitude is still insignificant ($t<t_{2}$). Here we determine what occurs if the encounter takes place well before $t_{1}$. For this purpose, we performed B-1-500 run, in which the fully scaled satellite approaches the host galaxy 500 Myr earlier than in the B-1 run.

In the surface density maps (Fig. 15), we see disc perturbations due to the encounter immediately after $t=0.35$ Gyr. While there is no visible random perturbation in the disc (see upper left panel), the only perturbation that exists is the one due to the satellite and subsequent delayed response in the central region ($R<2$ kpc), visible only after $t=0.4$ Gyr. Then it disappears at 0.6 Gyr, and a second perturbation appears after a short gap. Once the first random wave has reached the centre at $t\simeq 0.75$ Gyr, the large scale spiral becomes tightly wound and ineffective in providing a delayed response in the centre. This follows from comparison of the amplitudes in the left panels for $R<2$ kpc.

The resulting bar amplitude calculated for disc particles within characteristic scale length $R_{\textrm{d}}$ is shown in Fig. 13 together with the bar amplitude of the B-0 run. It is seen that after a jump at $t\simeq 0.4$ Gyr due to the delayed response waves, the central perturbation gradually fades away until the random wave forms a bar perturbation at $t\simeq 0.75$ Gyr – seen already in the previous runs. Afterwards, the B-1-500 curve continues almost in pace with the B-1 curve.

7 Comparison with observations

In this section we would like to discuss how our results compare with observations. Throughout this paper we have argued that the Aq-D2 is believed to be a fair representative of a Milky Way-like host. Fig. 3 in Moetazedian & Just (2016) is very helpful for comparing the distribution of satellites from the Aq-D2 run (red data points) with observations of our Milky Way. According to this figure, there resides a handful of massive satellites with $M_{\textrm{tid}}>10^{9}M_{\odot}$ among the seven simulations with pericentre distances within 25 kpc from the centre. The Milky Way’s Sagittarius dwarf galaxy (Sgr) has an estimated current total mass $M_{\textrm{Sgr}}\sim 5\times 10^{8}M_{\odot}$ and a pericentre of $\sim 20$ kpc (Law, Johnston & Majewski, 2005). If we assume Sgr to be on its second passage, then it has probably lost $\sim$50% of its mass during the initial crossing. In the Aq-D2, an analogue to Sgr can be found at $\sim 20$ kpc. In addition, more massive Milky Way satellites such as SMC with a mass $\sim 6.5\times 10^{9}M_{\odot}$ (Bekki & Stanimirović, 2009) and galactocentric distance of 60 kpc (McConnachie, 2012) have a corresponding analogue in the Aq-D2 realisation.

As previously shown by Eq. 2.10, the impact of a satellite with a given tidal mass and a pericentric distance $b$ scales as $\sim M_{\textrm{tid}}/b^{2}$. For a Sgr or SMC-like satellite we expect smaller impact as their passage occurs at a further distance to the centre (50-60 kpc) and hence, we fall below $10^{-4}$ (Fig. 2).

The distribution of satellites from both N-body cosmological simulations and observational studies are dominated by lower mass objects. In the Aquarius satellites distribution there exists a fair number of such low-mass satellites with $M_{\textrm{tid}}<10^{9}M_{\odot}$ and pericentres within the inner 10 kpc of the disc. Observations of the Local Group have shown the existence of so-called ultra-faint dwarf galaxies with velocity dispersions of the order of 10 kms*-1* (Simon & Geha, 2007). Recent observations as part of the Dark Energy Survey (DES) have increased the total number of detected ultra-faint dwarfs to 17 (Drlica-Wagner et al., 2015). This acts in favour of narrowing down the gap between cosmological model predictions and observations at the lower mass end of the galaxy mass function, that is, “the missing satellite problem”. The DES year-two quick release (Y2Q1) predicts a total of $\sim$100 ultra-faint dwarfs in the full sky cover. These objects are highly dark matter dominated with a very low stellar budget. In order to estimate the DM mass content of the satellites, their velocity dispersion needs to be measured using spectroscopy; this has only been done for a few satellites (e.g. Martin et al. 2016; Simon et al. 2016). These galaxies would be subject to strong tidal stripping in case of very close encounters with the disc, meaning it is very difficult to observe any remnants. Therefore, our current observational sample of ultra-faint dwarfs is limited and fainter magnitudes need to be reached.

There exists an inconsistency in current models of the Milky Way bar, where a broad range of values are discussed for bar properties such as age, pattern speed etc. (e.g. Wegg, Gerhard & Portail 2015; Monari et al. 2017). Therefore we are unable to constrain our models, making it impossible to determine the age of the Milky Way’s bar.

8 Conclusions

In this paper we analyse how the bar formation time can be affected by dark matter satellite encounters. Our Milky Way-like galaxy model is studied mainly using N-body simulations. All components, including an exponential stellar disc, a NFW halo, and a Sérsic cuspy bulge are tailored to the Milky Way and are live, that is, represented by particles that interact according to Newton’s law of gravity.

The adopted parameters of the host galaxy provide a disc mildly sensitive to bar instability with a characteristic e-folding time (inverse growth rate) $\tau_{\textrm{e}}=340$ Myr. This time, as well as the bar pattern speed, $\Omega_{\textrm{p}}=48...49$ km s*-1kpc-1* were obtained from our simulations and agree well with global mode calculations in the framework of linear perturbation theory. Such a large e-folding time is primarily due to the relatively high radial velocity dispersion (35 km s*-1* in the Solar neighbourhood), and a weak central cusp with $\rho\propto r^{-0.5}$. The cusp does not prohibit bar formation, as long as a disc with finite thickness rather than a razor-thin disc is considered (Polyachenko et al., 2016a).

Bar formation in the isolated galaxy occurs in the presence of random swing amplified wave packets appearing in a range of radii with Toomre parameters $X<3$ and $Q<3$ (between 2 and 4 kpc in our case). These packets interfere with a bar-like perturbation being a density wave in the inner part, which grows on its own due to bar instability. The interference leads randomly to demolition or enhancement of the bar. The latter is seen as a jump of the bar amplitude at some arbitrary moment in time ($t_{1}$ in Fig. 5). However, the exponential growth occurs later after additional random wave packets reinforce the inner bar, thus the moment $t_{2}$ at which the bar starts to grow exponentially is also a random quantity (see also Polyachenko, 2016).

The properties and orbital parameters of the satellite galaxies are derived from the Aquarius D-2 cosmological simulation. All satellites with tidal mass $M_{\textrm{tid}}$ above $10^{8}M_{\odot}$ and pericentric distances to the host galaxy below 25 kpc were selected. Our set of runs with satellites consists of four runs with one satellite, two runs with four satellites (with tidal mass M${}_{\textrm{tid}}\geq 10^{9}M_{\odot}$) and one run with all 23 satellites. Each satellite was represented by 50,000 particles. In all cases, the pattern speed and the e-folding time of the bar in the host galaxy remain unchanged. However, we noted that from run to run, the bar amplitudes grow asynchronously. In the earliest case, the bar forms 300 Myr earlier than in the reference isolated run without satellites. In the latest run, the bar appeared with a 600 Myr delay. Most of the advancement (or delay in the other two cases) is caused by the primary satellite whose pericentre is the closest to the galactic centre. The effect from other massive satellites is appreciable, yet turns out to be less significant. The lighter satellites play practically no role.

Satellites add complexity to the physical picture of bar formation. One satellite on the same orbit can give both advancement and delay, simply by passing in a slightly different time. Thus, orientation of the bar is important. However, it is difficult to specify ‘a rule of thumb’ for the sign of the effect, since the satellite angular speed is at least two times higher than the bar pattern speed, so a range of various bar phases is seen in all cases.

The main features of the satellite impact on the bar formation can be understood again through the wave interactions. After the satellite passage, a tidal wave is induced across the disc (Fig. 14), which then evolves in accordance with the WKB relation (5.21). It winds up and propagates to the centre where it interferes with the bar, as follows from Fig. 11. While it remains relatively open, it produces a delayed response in the inner disc. This mechanism is similar to swing amplification observed in simulations of the stable Mestel disc that turns out to be surprisingly vital under the influence of minor dipole perturbations (T81). This interference, in some cases, leads to bar formation delay (B-1m, B-1-110), or advancement (B-1), depending on enhancement or attenuation of the central structure. Keeping the phase of the tidal wave fixed, the interference pattern depends on the bar phase, which is illustrated by comparison of B-1 and B-1-110 runs. This basic mechanism is contaminated by additional random waves, which appear throughout the disc, also producing delayed responses in the inner disc, regardless of any encounter.

The effect is possible only if the encounter occurs after a seed bar is formed (after the amplitude jump at $t=t_{1}$) and before the exponential growth phase, $t\gtrsim t_{2}$. If the encounter occurs well before $t_{1}$, the tidal wave decays with no effect. This is proved by our B-1-500 run, in which the tidal wave occurs at $t\approx 400$ Myr, as follows from Fig. 15. The bar formation time then coincides with one in the isolated run B-0. During the period of vulnerability, $t_{1}\lesssim t\lesssim t_{2}$, the tidal wave amplitude is comparable to the amplitude of the seed bar, so the consequences of the interference are most evident. After $t_{2}$, the bar becomes strong so that it cannot be affected by relatively weak tidal waves.

It was known that the bar formation time depends on the instability properties of the galactic disc, which is determined by many parameters including the disc mass, velocity dispersion, dark matter halo, and a central matter concentration. This leads to uncertainties in predicting the bar age from simulations. The discussed phenomenon adds a scatter in the age of the bars. Polyachenko (2016) showed that less unstable models have longer time between the first jump and the start of the exponential growth. This implies that at higher redshifts $z>0.5$ satellites have more chances to interfere with the process of bar formation an thus the effect should be stronger. This analysis supports the hypothesis that all observed barred galaxies should be bar unstable in order to develop long-lasting bars.

Acknowledgements.

We would like to thank the anonymous referee for their constructive comments and remarks which have greatly improved this paper. The main production runs were done on the MilkyWay supercomputer, funded by the Deutsche Forschungsgemeinschaft (DFG) through the Collaborative Research Centre (SFB 881) “The Milky Way System” (subproject Z2), hosted and co-funded by the Jülich Supercomputing Center (JSC). The special GPU accelerated supercomputer laohu funded by NAOC/CAS and through the “Qianren” special foreign experts program of China (Silk Road Project), has been used for some of code development. We also used a smaller GPU cluster kepler for data analysis, funded under the grants I/80 041-043 and I/81 396 of the Volkswagen Foundation. This work was supported by the Sonderforschungsbereich SFB 881 “The Milky Way System” (subproject A2 and A6) of the German Research Foundation (DFG). E.P. acknowledges financial support by the Russian Basic Research Foundation, grants 15-52-12387, 16-02-00649, and by the Basic Research Program OFN-15 ‘The active processes in galactic and extragalactic objects’ of Department of Physical Sciences of RAS. P.B. acknowledges the special support by the NASU under the Main Astronomical Observatory GRID/GPU “golowood” computing cluster project.