© ESO, 2017

1. Introduction

Structure formation in the Universe governed by the Λ cold dark matter (Λ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 ≳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 et al. 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 et al. 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 et al. 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_sat/M_galaxy = 1 or 3 were able to only produce transient bars. Also, Walker et al. (1996) showed that a satellite with M_sat/M_disc = 0.1 on a circular prograde orbit with a pericentre distance of 6R_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_sat/M_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_d = 3.53 × 10¹⁰M_⊙ and R_d = 2.82 kpc, respectively. Also, the satellites were on eccentric orbits and possessed a mass range 0.21−0.57 M_d with pericentre distances of 0.5−6.2 R_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⁸M_⊙ extracted from cosmological simulations likely to host Milky Way-like galaxies (Springel et al. 2008; Moetazedian & Just 2016).

In Sect. 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 Sect. 4. The relevance of the azimuthal phase between a weak bar and incoming satellites towards the bar formation epoch is investigated in Sect. 5.1. Also, the results from the run with an early satellite encounter are presented in Sect. 6. Section 7 contains comparisons with observations. We finalise the paper with a summary and discussion of the results in Sect. 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_max, the radial distance of this velocity r_vmax and the tidal mass M_tid for every substructure (satellite) within the simulation box. The parent host halo has an enclosed mass, M₂₀₀ = 1.774 × 10¹²M_⊙, which corresponds to the mass within a sphere with radius r₂₀₀ = 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 (ρ_crit). Also the radius, where the halo encloses the mass with mean density 50 times ρ_crit is denoted as r₅₀ 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 Table 1; N_sub is the total number of satellites with M_tid ≥ 10⁶M_⊙ and f_sub (<r₅₀) is the fraction of these objects within r₅₀. 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_cross (<r₅₀). The last four rows show the number of crossed satellites with masses 1, 3, 5 and 10 × 10⁸M_⊙. We argue in Moetazedian & Just (2016) that only satellites with potentially contribute towards the heating of the Galactic disc. The Aq-D2 simulation has 23 satellites with out of which four have masses larger than 10⁹M_⊙.

The host halo, in which these satellites would be inserted, has a slightly different enclosed mass (M₂₀₀). 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_{200,Aq - D2}/M₂₀₀ (=1.37 in our case),

$\begin{array}{ccc}& & \\ & & \\ & & \end{array}$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) ${\mathit{\rho }}_{\mathrm{NFW}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}\frac{{\mathit{\rho }}_{\mathrm{s}}}{\mathrm{(}\mathit{r}\mathit{/}{\mathit{r}}_{\mathrm{s}}\mathrm{)}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{r}\mathit{/}{\mathit{r}}_{\mathrm{s}}{\mathrm{)}}^{\mathrm{2}}}$(4)and${\mathit{\delta }}_{\mathit{c}}\mathrm{=}\frac{{\mathit{\rho }}_{\mathrm{s}}}{{\mathit{\rho }}_{\mathrm{crit}}}\mathrm{=}\frac{\mathrm{200}}{\mathrm{3}}\frac{{\mathit{c}}^{\mathrm{3}}}{\ln \mathrm{(}\mathrm{1}\mathrm{+}\mathit{c}\mathrm{)}\mathrm{-}\mathit{c}\mathit{/}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{c}\mathrm{)}}\mathrm{·}$(5)Here ρ_s and r_s are the scale density and the scale radius of the halo, respectively. The density contrast δ_c with respect to ρ_crit is calculated with the concentration c of the halo being the ratio of r₂₀₀/r_s.

As discussed in Sect. 2.3, all the DM satellites also follow a NFW profile, which is tidally truncated at r_tid according to the satellite mass M_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 Table 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_d = 2.9 kpc and truncation radius 15 kpc. The radial velocity dispersion σ_R is approximately exponential with central value σ_R0 = 140 km s^-1 and radial scale length R_σ = 2R_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${\mathit{\rho }}_{\mathrm{b}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}{\mathit{\rho }}_{\mathrm{b}}{\left(\frac{\mathit{r}}{{\mathit{R}}_{\mathrm{e}}}\right)}^{\mathrm{-}\mathit{p}}{\mathrm{e}}^{\mathrm{-}\mathit{b}\mathrm{(}\mathit{r}\mathit{/}{\mathit{R}}_{\mathrm{e}}{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathit{n}}}\mathit{,}$(6)where r is the spherical radius, corresponding to a Sérsic surface brightness profile with effective radius R_e. The scale density ρ_b, can be replaced by the bulge velocity scale ${\mathit{\sigma }}_{\mathrm{b}}\mathrm{\equiv }{\mathrm{\left\{}\mathrm{4}\mathit{\pi Gn}{\mathit{b}}^{\mathit{n}\mathrm{(}\mathit{p}\mathrm{-}\mathrm{3}\mathrm{)}}\mathrm{\Gamma }\mathrm{\right[}\mathit{n}\mathrm{(}\mathrm{3}\mathrm{-}\mathit{p}\mathrm{)}\mathrm{\left]}{\mathit{R}}_{\mathrm{e}}^{\mathrm{2}}{\mathit{\rho }}_{\mathrm{b}}\mathrm{\right\}}}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathit{.}$(7)Here, ${\mathit{\sigma }}_{\mathrm{b}}^{\mathrm{2}}$ corresponds to the depth of the gravitational potential associated with the bulge. In addition, we require n = 1.11788 leading to p ≃ 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₂₀₀. The bulge component dominates the inner rotation curve at R ≳ 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 ≲ 1.5 kpc, while the halo takes over at R> 10 kpc. At R ≈ 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.

Fig. 1 Initial profiles for the basic model: a) the total circular velocity and its components due to disc, bulge, and halo; b) angular velocity Ω(R) and curves Ω(R) ± κ(R) / 2; c) Toomre Q initially set by GalactICS and the one actually obtained in the simulations.

Panel (b) presents the angular velocity profile Ω(R) in the equatorial plane. In addition, positions of the inner (ILR) and outer Lindblad resonances (OLR) are shown as Ω ± κ/ 2. For a given pattern speed Ω_p, the ILR is calculated using${\mathrm{\Omega }}_{\mathrm{p}}\mathrm{=}{\mathrm{\Omega }}_{\mathrm{pr}}\mathrm{\left(}\mathit{R}\mathrm{\right)}\mathrm{\equiv }\mathrm{\Omega }\mathrm{\left(}\mathit{R}\mathrm{\right)}\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}\mathit{\kappa }\mathrm{\left(}\mathit{R}\mathrm{\right)}\mathit{.}$(8)The quantity, Ω_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 → 0 with R^{− α/ 2} and α ≈ 0.5.

The lower panel (c) of Fig. 1 demonstrates the Toomre Q profile $\mathit{Q}\mathrm{=}\frac{\mathit{\kappa }{\mathit{\sigma }}_{\mathit{R}}}{\mathrm{3.36}\mathit{G}{\mathrm{\Sigma }}_{\mathrm{d}}}\mathit{,}$(9)where Σ_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_min = 1.8 resides at R = 5.9 kpc. This minimum rises to Q_min ≈ 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 ρ_s and r_s in Eq. (4)are determined by r_vmax using r_s = 0.4623 r_vmax and V_max using the enclosed mass inside r_vmax. Next we derive the tidal radius r_tid as cut-off radius from the cumulative mass profile reaching the satellite mass M_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 → ∞; therefore, we use the satellites’ r_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_c(R_d): ${\left(\frac{\mathrm{\Delta }\mathit{v}}{{\mathit{v}}_{\mathrm{c}}}\right)}_{{\mathit{R}}_{\mathrm{d}}}\mathrm{~}\frac{\mathrm{2}\mathit{G}{\mathit{M}}_{\mathrm{tid}}}{{\mathit{b}}^{\mathrm{2}}\mathit{V}}\frac{{\mathit{R}}_{\mathrm{d}}}{{\mathit{v}}_{\mathrm{c}}}\mathrm{~}\frac{\mathit{G}{\mathit{M}}_{\mathrm{tid}}{\mathit{R}}_{\mathrm{d}}}{{\mathit{b}}^{\mathrm{2}}{\mathit{v}}_{\mathrm{c}}^{\mathrm{2}}}\mathrm{·}$(10)Here b is the pericentre distance and we assumed that velocities at the time of encounters V are similar and approximated as V ≃ 2v_c. Figure 2 shows the reverse cumulative histogram of the number of encounters with the tidal impact parameter (10) greater than a given one, for all encounters with M_tid> 10⁸M_⊙ 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 ~ 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 ~ 2, compared to the most massive satellite.

Fig. 2 The reverse cumulative histograms of the number of encounters with the tidal impact parameter (10) greater than the given one, for all seven runs of Moetazedian & Just (2016). We consider the encounters with Mtid> 108M⊙ (black) and Mtid> 109M⊙ (blue) before applying the rescaling. The red arrow marks the position of the Aq-D2 encounter used in simulations with one satellite and it is the 5th encounter from the right hand-side of the upper cumulative line. The satellite with the highest impact (~ 0.4) has a mass 1.1 × 108M⊙ and comes as close as 0.25 kpc to the centre.

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_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 Table 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,${Ax}\mathrm{=}\mathit{\omega }{x}\mathit{,}$(11)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 Ω_b ≈ 48.34 km s^-1 kpc^-1 and for the growth rate we obtained ω_I ≈ 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-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 θ = 0.5, accompanied by individual gravitational softening of 10 pc. bonsai2 employs the leap-frog integration scheme with a fixed time–step Δ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.

Table 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. (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_tid = 2.45 × 10⁹M_⊙. For runs marked as “B-4”, the four most massive Aq-D2 crossed satellites with M_tid> 10⁹M_⊙ were employed.

Fig. 3 The orbit of the primary satellite from our five runs as a function of time.

The orbital positions and velocities of the satellites are determined by their density centres. Figure 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.

Fig. 4 Bar patterns in the B-0 run oriented along the x-axis at different stages of bar evolution. The curves are isolines of the surface density evenly spaced in log scale (ten levels for every factor of 10). Each frame size is 32 × 32 kpc.

The bar is a wave, and its growth can be investigated using relative amplitudes of perturbed quantities, such as the surface density A₂/A₀,${\mathit{A}}_{\mathit{m}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{=}\sum _{\mathit{j}}{\mathit{\mu }}_{\mathit{j}}{\mathrm{e}}^{\mathrm{-}\mathrm{i}\mathit{m}{\mathit{\theta }}_{\mathit{j}}}\mathit{,}$(12)(here μ_j and θ_j are mass and polar angle of star j; j spans particles within some fixed radius, for example, the radial scale length R_d), or velocity components U_k,m. The latter are defined as the perturbed velocity amplitudes,${\mathit{V}}_{\mathit{k,m}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{{\mathit{A}}_{\mathrm{0}}}\sum _{\mathit{j}}{\mathit{\mu }}_{\mathit{j}}\begin{array}{c}\mathrm{˜}\\ {\mathit{v}}_{\mathit{k,j}}\end{array}{\mathrm{e}}^{\mathrm{-}\mathrm{i}\mathit{m}{\mathit{\theta }}_{\mathit{j}}}\mathit{,}$(13)normalized to the average velocity dispersion σ_k(t) calculated by${\mathit{\sigma }}_{\mathit{k}}^{\mathrm{2}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{{\mathit{A}}_{\mathrm{0}}}\sum _{\mathit{j}}{\mathit{\mu }}_{\mathit{j}}\stackrel{}{{{\mathit{v}}_{\mathit{k,j}}^{\mathrm{2}}}_{\mathrm{˜}}}\mathit{.}$(14)The cylindrical components R, θ, and z correspond to the quantity k, and $\begin{array}{c}{\mathrm{˜}}_{}\\ {\mathit{v}}_{\mathit{k,j}}\end{array}$ is the residual velocity component of the particle j.

Panel (a) in Fig. 5 shows the relative bisymmetric amplitudes A₂/A₀, 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 Table 3, t₁ denotes the jump, t₂ corresponds to the start of the exponential growth and t₃ 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, $\mathit{B}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{=}\mathrm{1}\mathrm{-}{\mathit{I}}_{\mathit{yy}}\mathit{/}{\mathit{I}}_{\mathit{xx}}\mathit{,}$(15)give estimates for the growth rates, which appear to be very close to one another (and ≃ 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 Ω_p determined by an average centred finite difference over two time-intervals of ± 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₁. Just after the jump, the centred finite difference shows much less scatter, and the filtered Ω_p varies only slowly until t = 1.9 Gyr. Due to a bar slowdown, the pattern speed gradually decreases to 28.6 km s^-1 kpc^-1 at t = 4 Gyr.

Fig. 5 Bar formation in the isolated model. a) Bar amplitudes A2/A0 and URθ,2, and the bar strength B(t) demonstrate lags, jumps, and exponential growths with a rate ωI ≈ 3 Gyr-1 (the slope is shown by the black dashed curve). b) The pattern speed of the bar Ωp determined as an average centred finite difference of the bar phase (dots), and the spline smoothed Ωp (red curve).

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_b, as a radius at which the ellipticity ε(r) ≡ 1−b_e/a_e (a_e and b_e are the major and minor semi-axes of ellipses) declines by ~ 15% from its maximal value. The obtained radii and ellipticities are given in Table 4. For t = 1.0 Gyr, the bar radius is approximately 1.1 kpc, and ε ≃ 0.13. At the end of the exponential growth, t ≈ 2.5, the bar radius is 4.1 kpc, with an ellipticity of ε ≃ 0.65.

With the determined Ω_p ≈ 49 km s^-1 kpc^-1 at T = 1 Gyr, the corotation resonance occurs at R_C ≈ 4.4 kpc. The ILR and OLR are located at 0.7 and 7.9 kpc, respectively. As the bar grows between t₂ and t₃, the ratio ℛ ≡ R_C/R_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₁, 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₁. However, the amplitude growth now is irregular and slower than exponential, until the wave amplitude A₂ reaches 1 or 2 per cent of the axisymmetric background A₀ (t = t₂). The exponential growth with a growth rate ω_I ≈ 3 Gyr^-1 lasts until t = t₃, 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₂/A₀ of these models with the one calculated for B-0 (solid red). The jumps in all curves occur at the same time, t₁ = 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 Δ_B is given in Table 3.

Fig. 6 Bar formation in models with satellites: a) bar amplitudes A2/A0 are given for B-0, B-1, B-1m, B-4, B-4m and B-23 run. The black dashed line shows the slope corresponding to a growth rate ωI = 3 Gyr-1; b) the spline smoothed pattern speeds for the same models.

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 Δ_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 Sect. 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 ≈ 400 km s^-1; for a galaxy of size d ~ 20 kpc the encounter lasts d/V ≈ 50 Myr, which is a fraction of the rotation period of the bar, ≈ 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: $\mathrm{d}{v}\mathrm{\left(}\mathit{R,\theta }\mathrm{\right)}\mathrm{=}\mathrm{-}\mathit{G}\frac{{\mathit{M}}_{\mathrm{s}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{(}{R}\mathrm{-}{{r}}_{\mathrm{s}}\mathrm{)}}{\mathrm{|}{{r}}_{\mathrm{s}}\mathrm{-}{R}{\mathrm{|}}^{\mathrm{3}}}\mathrm{d}\mathit{t,} \mathrm{d}\mathit{\theta }\mathrm{=}\mathrm{\Omega }\mathrm{\left(}\mathit{R}\mathrm{\right)}\mathrm{d}\mathit{t,}$(16)where M_s(t) and r_s(t) are known functions for the satellite mass and its distance from the disc centre, 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_θ,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.

Fig. 7 Time sequence of phase (upper panels) and amplitudes (lower panels) of the m = 2 perturbations of radial (red triangles) and azimuthal (blue circles) velocities induced by the primary satellite (pink diamonds, visible for T = 0.84−0.86 Gyr) in the cold fluid disc.

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 ∝ e^i(mθ−ωt), the perturbed surface density Σ_da and average velocity components (v_θa,v_Ra) obey the following relations (Binney & Tremaine 2008):

$\begin{array}{ccc}& & \\ & & \\ & & \end{array}$where B ≡ −κ²/ 4Ω is the Oort constant; Σ₀ is the disc surface density and Φ_a is the potential perturbation; $\mathrm{\Delta }\mathrm{\equiv }{\mathit{\kappa }}^{\mathrm{2}}\mathrm{-}\mathrm{(}\mathit{\omega }\mathrm{-}\mathit{m}\mathrm{\Omega }{\mathrm{)}}^{\mathrm{2}}$(20)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 ${\mathrm{\Phi }}_{\mathit{a}}\mathrm{\left(}\mathit{R}\mathrm{\right)}\mathrm{=}\mathrm{\left|}{\mathrm{\Phi }}_{\mathit{a}}\mathrm{\right|}{\mathrm{e}}^{\mathrm{i}\mathit{F}\mathrm{\left(}\mathit{R}\mathrm{\right)}}\mathrm{=}\mathrm{\left|}{\mathrm{\Phi }}_{\mathit{a}}\mathrm{\right|}{\mathrm{e}}^{\mathrm{i}{\mathrm{\int }}^{\mathit{R}}\mathit{k}\hspace{0.17em}\mathrm{d}\mathit{R}}\mathit{,}$(21)with phase function F(R), where k = dF(R) / dR and | kR | ≫ 1. From Eq. (19)one has $\mathit{k}{\mathit{v}}_{\mathit{Ra}}{\mathrm{\Sigma }}_{\mathrm{0}}\mathrm{=}\mathrm{(}\mathit{\omega }\mathrm{-}\mathit{m}\mathrm{\Omega }\mathrm{)}{\mathrm{\Sigma }}_{\mathit{da}}\mathit{,}$(22)that is, surface density of the wave is in phase with the radial velocity outside corotation, and in anti-phase inside corotation: $\begin{array}{c}\\ {\mathit{F}}_{\mathit{R}}\mathrm{-}{\mathit{F}}_{\mathrm{\Sigma }}\mathrm{=}\mathit{\pi }& \mathrm{:}& \mathit{R}\mathit{<}{\mathit{R}}_{\mathit{C}}\mathit{,}\\ {\mathit{F}}_{\mathit{R}}\mathrm{-}{\mathit{F}}_{\mathrm{\Sigma }}\mathrm{=}\mathrm{0}& \mathrm{:}& \mathit{R}\mathit{>}{\mathit{R}}_{\mathit{C}}\mathit{.}\end{array}$(23)The equation for the azimuthal velocity (Eq. (18)) simplifies to ${\mathit{v}}_{\mathit{\theta a}}{\mathrm{\Sigma }}_{\mathrm{0}}\mathrm{=}\frac{\mathrm{-}\mathrm{2}\mathrm{i}\mathit{B}}{\mathrm{\Delta }}\mathit{k}{\mathrm{\Phi }}_{\mathit{a}}\mathit{,}$(24)and taking into account that Φ_a = −2πGΣ_da/ | k |, ${\mathit{F}}_{\mathit{\theta }}\mathrm{-}{\mathit{F}}_{\mathrm{\Sigma }}\mathrm{=}\mathrm{-}\mathit{\pi }\mathit{/}\mathrm{2}\mathit{,}$(25)meaning that the surface density lags behind the azimuthal velocity by π/ 4 for an m = 2 perturbation.

Fig. 8 Loci of maxima of m = 2 perturbations (velocity and surface density) in the B-0 run (no satellites) at the time of the encounter in B-1 (the pink diamond shows the position of the primary satellite in B-1 for visualisation).

In the bar region we assume that the radial derivative is negligible compared to the other term in square brackets in Eqs. (17), (18)). For the radial velocity, ${\mathit{v}}_{\mathit{Ra}}\mathrm{\left(}\mathit{R}\mathrm{\right)}\mathrm{=}\mathrm{-}\frac{\mathrm{2}\mathrm{i}\mathit{m}}{\mathrm{\Delta }}\frac{\mathrm{\Omega }}{\mathit{R}}{\mathrm{\Phi }}_{\mathit{a}}\mathit{,}$(26)that is, ${\mathit{F}}_{\mathit{R}}\mathrm{-}{\mathit{F}}_{\mathrm{\Sigma }}\mathrm{=}\mathit{\pi }\mathit{/}\mathrm{2}\mathit{,}$(27)meaning that radial velocity lags behind the surface density by π/ 4 (m = 2). For the azimuthal velocity, ${\mathit{v}}_{\mathit{\theta a}}\mathrm{\left(}\mathit{R}\mathrm{\right)}\mathrm{=}\mathrm{-}\frac{\mathit{m}}{\mathrm{\Delta }}\frac{\mathrm{(}\mathit{\omega }\mathrm{-}\mathit{m}\mathrm{\Omega }\mathrm{)}}{\mathit{R}}{\mathrm{\Phi }}_{\mathit{a}}\mathit{,}$(28)that is, $\begin{array}{c}\\ {\mathit{F}}_{\mathit{\theta }}\mathrm{-}{\mathit{F}}_{\mathrm{\Sigma }}\mathrm{=}\mathit{\pi }& \mathrm{:}& \mathit{R}\mathit{<}{\mathit{R}}_{\mathit{C}}\mathit{,}\\ {\mathit{F}}_{\mathit{\theta }}\mathrm{-}{\mathit{F}}_{\mathrm{\Sigma }}\mathrm{=}\mathrm{0}& \mathrm{:}& \mathit{R}\mathit{>}{\mathit{R}}_{\mathit{C}}\mathit{.}\end{array}$(29)

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₂(R,t) /A₀(R,t), ${\mathit{A}}_{\mathit{m}}\mathrm{\left(}\mathit{R,t}\mathrm{\right)}\mathrm{=}{\sum _{\mathit{j}}}^{\mathrm{\prime }}{\mathit{\mu }}_{\mathit{j}}{\mathrm{e}}^{\mathrm{-}\mathrm{i}\mathit{m}{\mathit{\theta }}_{\mathit{j}}}\mathit{,}$(30)and velocity components, V_k,2(R,t) /A₀(R,t), ${\mathit{V}}_{\mathit{k,m}}\mathrm{\left(}\mathit{R,t}\mathrm{\right)}\mathrm{=}{\sum _{\mathit{j}}}^{\mathrm{\prime }}{\mathit{\mu }}_{\mathit{j}}\begin{array}{c}\mathrm{˜}\\ {\mathit{v}}_{\mathit{k,j}}\end{array}{\mathrm{e}}^{\mathrm{-}\mathrm{i}\mathit{m}{\mathit{\theta }}_{\mathit{j}}}\mathit{.}$(31)In both expressions, the summations are taken over disc particles within a ring ΔR near R.

Figure 8 shows the loci of maxima θ(R) = −F/m for the B-0 run without satellites. Note that the maxima of V_θ,2 are shifted by π/ 2 from the surface density maxima, and the maxima of V_R,2 lag by π/ 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.

Figure 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 ~ 110 Myr). So in this case the bar is out of phase when the satellite approaches the pericentre.

Fig. 9 Same as in Fig. 8, but at the time of the encounter in B-1m.

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₂/A₀ | 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).

Fig. 10 Bar amplitudes | A2/A0 | for the B-0, B-1, B-1m runs (selected zoom of Fig. 6).

Fig. 11 Surface density maps log | A2/A0 | (left panels, red = high, blue = low amplitude); Re A2/A0 (right panels, red = positive, blue = negative amplitude at the x-axis) for B-0, B-1, B-1m, and B-1-110 runs (top to bottom), reflecting the interference between the bar-like perturbation in the centre, random perturbations, and tidally induced waves. The time period is 0.6 <t< 1.8 Gyr.

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~1...2 R_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 Sect. 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, Ω(R). This is a manifestation of “a strong cooperative effect” called swing amplification. For this mechanism to be efficient, the parameter$\mathit{X}\mathrm{\equiv }\frac{{\mathit{\kappa }}^{\mathrm{2}}\mathit{R}}{\mathrm{2}\mathit{m\pi G}{\mathrm{\Sigma }}_{\mathrm{0}}}$(32)should not exceed 3, which is the case in our model in the range 0.6 <R< 5.5 kpc, with a minimum X ≃ 2 at R ≃ 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 Re A₂/A₀, 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 ~ 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_i = aR + b. Since frequency difference between the bar and the spiral response Ω_b−Ω_s(R) = 2π/T_i, one can calculate the local pattern speed Ω_s(R) for wave patterns due to the encounter. This is shown by red circles in Fig. 12, along with Ω and Ω−κ/ 2 curves.

Fig. 12 Local pattern speeds Ωs of the wave patterns (red circles) compared to main frequencies of the disc (Ω is the angular velocity, κ is the epicyclic frequency, Ω∗ is given by Eq. (36)).

Local pattern speeds Ω_s of the wave patterns can be explained using the simplest fluid WKB dispersion relation. Each wave pattern evolves according to the relation,$\mathit{\omega }\mathrm{=}\mathit{\omega }\mathrm{\left(}\mathit{k,R}\mathrm{\right)}\mathit{,}$(33)meaning that R = R(k,ω) is a function of k; ω ≡ mΩ_s is a parameter. For very open spiral², k is close to zero, so Ω_s = Ω(R)−κ(R) /m, and R is close to the ILR, R_ILR. As the local wave pattern evolves, (ω−mΩ)² 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: $\mathrm{(}\mathit{\omega }\mathrm{-}\mathit{m}\mathrm{\Omega }\mathrm{(}\mathit{R}\mathrm{)}{\mathrm{)}}^{\mathrm{2}}\mathrm{=}{\mathit{\kappa }}^{\mathrm{2}}\mathrm{(}\mathit{R}\mathrm{)}\mathrm{-}\mathrm{2}\mathit{\pi G}{\mathrm{\Sigma }}_{\mathrm{0}}\mathrm{|}\mathit{k}\mathrm{|}\mathrm{+}{\mathit{k}}^{\mathrm{2}}{\mathit{c}}^{\mathrm{2}}\mathit{,}$(34)where k is the wavenumber and c is a sound speed, which we use interchangeably with the radial velocity dispersion. Thus ${\mathit{k}}_{\mathrm{\ast }}\mathrm{=}\frac{\mathit{\pi G}{\mathrm{\Sigma }}_{\mathrm{0}}}{{\mathit{c}}^{\mathrm{2}}}\mathit{,}$(35)and substitution to the WKB relation gives ${\mathrm{\Omega }}_{\mathrm{\ast }}\mathrm{\left(}\mathit{R}\mathrm{\right)}\mathrm{=}\mathrm{\Omega }\mathrm{\left(}\mathit{R}\mathrm{\right)}\mathrm{-}\frac{\mathit{\kappa }\mathrm{\left(}\mathit{R}\mathrm{\right)}}{\mathit{m}}\mathrm{(}\mathrm{1}\mathrm{-}\mathit{Q}\mathrm{(}\mathit{R}{\mathrm{)}}^{-2}{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathit{.}$(36)In the last expression, Q denotes Toomre stability parameter for fluid discs, $\mathit{Q}\mathrm{\left(}\mathit{R}\mathrm{\right)}\mathrm{=}\frac{\mathit{\kappa }\mathrm{\left(}\mathit{R}\mathrm{\right)}\mathit{c}\mathrm{\left(}\mathit{R}\mathrm{\right)}}{\mathit{\pi G}{\mathrm{\Sigma }}_{\mathrm{0}}\mathrm{\left(}\mathit{R}\mathrm{\right)}}\mathrm{·}$(37)As is seen from Fig. 12, Ω_∗ and Ω_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₁), but its amplitude is still insignificant (t<t₂). Here we determine what occurs if the encounter takes place well before t₁. 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 ≃ 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_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 ≃ 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 ≃ 0.75 Gyr – seen already in the previous runs. Afterwards, the B-1-500 curve continues almost in pace with the B-1 curve.

Fig. 13 Bar amplitudes for B-0 and B-1-500 runs.

Fig. 14 Onset of the large-scale bar-like perturbation from the satellite in the B-1 run (time from top to bottom, Σd, VR and Vθfrom left to right).

Fig. 15 As in Fig. 11 for B-0 and B-1-500 runs for the period 0.2 < t < 1.4 Gyr.

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. Figure 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_tid > 10⁹M_⊙ 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_Sgr ~ 5 × 10⁸M_⊙ and a pericentre of ~ 20 kpc (Law et al. 2005). If we assume Sgr to be on its second passage, then it has probably lost ~50% of its mass during the initial crossing. In the Aq-D2, an analogue to Sgr can be found at ~ 20 kpc. In addition, more massive Milky Way satellites such as SMC with a mass ~ 6.5 × 10⁹M_⊙ (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. (10), the impact of a satellite with a given tidal mass and a pericentric distance b scales as ~ M_tid/b². 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_tid< 10⁹M_⊙ 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 ~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. 2017). 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 et al. 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) τ_e = 340 Myr. This time, as well as the bar pattern speed, Ω_p = 48...49 km s^-1 kpc^-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 ρ ∝ 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₁ in Fig. 5). However, the exponential growth occurs later after additional random wave packets reinforce the inner bar, thus the moment t₂ 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_tid above 10⁸M_⊙ 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_tid ≥ 10⁹M_⊙) 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 (36). 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₁) and before the exponential growth phase, t ≳ t₂. If the encounter occurs well before t₁, the tidal wave decays with no effect. This is proved by our B-1-500 run, in which the tidal wave occurs at t ≈ 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₁ ≲ t ≲ t₂, the tidal wave amplitude is comparable to the amplitude of the seed bar, so the consequences of the interference are most evident. After t₂, 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.

Acknowledgments

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.

References

All Tables

All Figures

	Fig. 1 Initial profiles for the basic model: a) the total circular velocity and its components due to disc, bulge, and halo; b) angular velocity Ω(R) and curves Ω(R) ± κ(R) / 2; c) Toomre Q initially set by GalactICS and the one actually obtained in the simulations.
In the text

Fig. 2 The reverse cumulative histograms of the number of encounters with the tidal impact parameter (10) greater than the given one, for all seven runs of Moetazedian & Just (2016). We consider the encounters with Mtid> 108M⊙ (black) and Mtid> 109M⊙ (blue) before applying the rescaling. The red arrow marks the position of the Aq-D2 encounter used in simulations with one satellite and it is the 5th encounter from the right hand-side of the upper cumulative line. The satellite with the highest impact (~ 0.4) has a mass 1.1 × 108M⊙ and comes as close as 0.25 kpc to the centre.

In the text

	Fig. 3 The orbit of the primary satellite from our five runs as a function of time.
In the text

	Fig. 4 Bar patterns in the B-0 run oriented along the x-axis at different stages of bar evolution. The curves are isolines of the surface density evenly spaced in log scale (ten levels for every factor of 10). Each frame size is 32 × 32 kpc.
In the text

	Fig. 5 Bar formation in the isolated model. a) Bar amplitudes A2/A0 and URθ,2, and the bar strength B(t) demonstrate lags, jumps, and exponential growths with a rate ωI ≈ 3 Gyr-1 (the slope is shown by the black dashed curve). b) The pattern speed of the bar Ωp determined as an average centred finite difference of the bar phase (dots), and the spline smoothed Ωp (red curve).
In the text

	Fig. 6 Bar formation in models with satellites: a) bar amplitudes A2/A0 are given for B-0, B-1, B-1m, B-4, B-4m and B-23 run. The black dashed line shows the slope corresponding to a growth rate ωI = 3 Gyr-1; b) the spline smoothed pattern speeds for the same models.
In the text

	Fig. 7 Time sequence of phase (upper panels) and amplitudes (lower panels) of the m = 2 perturbations of radial (red triangles) and azimuthal (blue circles) velocities induced by the primary satellite (pink diamonds, visible for T = 0.84−0.86 Gyr) in the cold fluid disc.
In the text

	Fig. 8 Loci of maxima of m = 2 perturbations (velocity and surface density) in the B-0 run (no satellites) at the time of the encounter in B-1 (the pink diamond shows the position of the primary satellite in B-1 for visualisation).
In the text

	Fig. 9 Same as in Fig. 8, but at the time of the encounter in B-1m.
In the text

	Fig. 10 Bar amplitudes | A2/A0 | for the B-0, B-1, B-1m runs (selected zoom of Fig. 6).
In the text

	Fig. 11 Surface density maps log | A2/A0 | (left panels, red = high, blue = low amplitude); Re A2/A0 (right panels, red = positive, blue = negative amplitude at the x-axis) for B-0, B-1, B-1m, and B-1-110 runs (top to bottom), reflecting the interference between the bar-like perturbation in the centre, random perturbations, and tidally induced waves. The time period is 0.6 <t< 1.8 Gyr.
In the text

	Fig. 12 Local pattern speeds Ωs of the wave patterns (red circles) compared to main frequencies of the disc (Ω is the angular velocity, κ is the epicyclic frequency, Ω∗ is given by Eq. (36)).
In the text

	Fig. 13 Bar amplitudes for B-0 and B-1-500 runs.
In the text

	Fig. 14 Onset of the large-scale bar-like perturbation from the satellite in the B-1 run (time from top to bottom, Σd, VR and Vθfrom left to right).
In the text

	Fig. 15 As in Fig. 11 for B-0 and B-1-500 runs for the period 0.2 < t < 1.4 Gyr.
In the text