Issue 
A&A
Volume 604, August 2017



Article Number  A75  
Number of page(s)  15  
Section  Galactic structure, stellar clusters and populations  
DOI  https://doi.org/10.1051/00046361/201630024  
Published online  09 August 2017 
Effects of galaxy–satellite interactions on bar formation
^{1} Zentrum für Astronomie der Universität Heidelberg, Astronomisches RechenInstitut, Mönchhofstr. 1214, 69120 Heidelberg, Germany
email: just@ari.uniheidelberg.de
^{2} Institute of Astronomy, Russian Academy of Sciences, 48 Pyatnitskya St., 119017 Moscow, Russia
^{3} Main Astronomical Observatory, National Academy of Sciences of Ukraine, MAO/NASU, 27 Akad. Zabolotnoho St. 03680 Kyiv, Ukraine
^{4} National Astronomical Observatories and Key Laboratory of Computational Astrophysics, Chinese Academy of Sciences, NAOC/CAS, 20A Datun Rd., Chaoyang District, Beijing 100012, PR China
Received: 7 November 2016
Accepted: 9 March 2017
Aims. We aim to show how encounters with lowmass satellite galaxies may alter the bar formation in a Milky Waylike disc galaxy.
Methods. We use highresolution Nbody simulations of a disc galaxy prone to mild bar instability. For realistic initial conditions of satellites, we take advantage of cosmological simulations of Milky Waylike 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 ~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_{⊙} 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 presentday host galaxy is expected to be even more significant at redshifts z ≳ 0.5.
Key words: galaxies: kinematics and dynamics / methods: numerical / Galaxy: disk
© 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ándezToledo et al. 2005; Lotz et al. 2010; Kannan et al. 2015). Current hydrodynamical cosmological simulations, which account for different small and largescale 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 nonaxisymmetric 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 Nbody 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 barunstable 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 zoomin cosmological simulations of Milky Waylike 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^{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 nonaxisymmetric (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 Waylike galaxy using highresolution Nbody 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_{⊙} extracted from cosmological simulations likely to host Milky Waylike 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. Nbody 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 DMonly 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 Waylike galaxies and which have not experienced a recent major merger. For the purpose of this study the level 2 of the AquariusD simulation, hereafter AqD2, was employed. We have shown in Moetazedian & Just (2016) that AqD2 is a fair representative of a typical Milky Waylike 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_{200} = 1.774 × 10^{12}M_{⊙}, which corresponds to the mass within a sphere with radius r_{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 (ρ_{crit}). Also the radius, where the halo encloses the mass with mean density 50 times ρ_{crit} is denoted as r_{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 Table 1; N_{sub} is the total number of satellites with M_{tid} ≥ 10^{6}M_{⊙} and f_{sub} (<r_{50}) is the fraction of these objects within r_{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_{cross} (<r_{50}). The last four rows show the number of crossed satellites with masses 1, 3, 5 and 10 × 10^{8}M_{⊙}. We argue in Moetazedian & Just (2016) that only satellites with potentially contribute towards the heating of the Galactic disc. The AqD2 simulation has 23 satellites with out of which four have masses larger than 10^{9}M_{⊙}.
Number statistics of AqD2 satellites.
The host halo, in which these satellites would be inserted, has a slightly different enclosed mass (M_{200}). Therefore, we needed to rescale the phasespace 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_{200} (=1.37 in our case),
The subscript “orig” represents the original nonscaled values. Early CDM Nbody simulations have shown that the density profile of DM (sub)haloes could be well fitted via the known NavarroFrankWhite (NFW) profile (Navarro et al. 1997) (4)and(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_{200}/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 threecomponent 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 BlandHawthorn & 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).
Characteristics of disk+bulge+halo for our host galaxy model.
The bulge component takes a density profile of the following form(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 (7)Here, 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_{200}. 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(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 (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 Nbody 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 cutoff 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 cutoff 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}): (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^{8}M_{⊙} among seven runs discussed in Moetazedian & Just (2016). The AqD2 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 AqD2 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 AqD2. 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 M_{tid}> 10^{8}M_{⊙} (black) and M_{tid}> 10^{9}M_{⊙} (blue) before applying the rescaling. The red arrow marks the position of the AqD2 encounter used in simulations with one satellite and it is the 5th encounter from the right handside of the upper cumulative line. The satellite with the highest impact (~ 0.4) has a mass 1.1 × 10^{8}M_{⊙} and comes as close as 0.25 kpc to the centre. 
The host halo has a different enclosed mass and concentration compared to the AqD2 main halo. This means the orbits of AqD2 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 AqD2 satellites using the factor f = 1.37, which corresponds to the ratio of AqD2 and the targeted host halo mass. Simulations with these initial conditions are called “fullyrescaled”. 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 razorthin 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 zdependence 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,(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. Nbody 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 Nbody simulations (Polyachenko et al. 2016a) we used three different codes (SUPERBOX–10, bonsai2, bergal0) and demonstrated the practical independence of the results on the chosen code. The simulations in this work are restricted to the bonsai2 treecode and bergal0^{1} (Zinchenko et al. 2015) routines as an auxiliary tool to analyse the snapshot data. The modified version of the recently developed Nbody 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 leapfrog 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 multimass 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 multimass and single mass runs varies from 10 to 100, thus the effective numerical resolution there is enhanced by this factor.
Summary of the runs.
The subscript “m” corresponds to the simulations, which are only massrescaled using Eq. (1); while the remaining simulations are fullyrescaled. In the case of simulations marked using “B1”, we use the primary satellite from AqD2 simulation with rescaled mass M_{tid} = 2.45 × 10^{9}M_{⊙}. For runs marked as “B4”, the four most massive AqD2 crossed satellites with M_{tid}> 10^{9}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 AqD2’s primary satellite from our five runs. There exists a 110 Myr difference between the first pericentre passage of B1 and B1m 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 B0 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 xaxis. 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 B0 run oriented along the xaxis 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_{2}/A_{0},(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,(13)normalized to the average velocity dispersion σ_{k}(t) calculated by(14)The cylindrical components R, θ, and z correspond to the quantity k, and 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 B0 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_{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, (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 timeintervals 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_{1}. 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 A_{2}/A_{0} and U_{Rθ,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 (MartinezValpuesta 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 semiaxes 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.
Parameters of the bar (bar radius, pattern speed, corotation radius, ellipticity) in B0 run at different moments in time.
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_{2} and t_{3}, 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 slowdown of the bar in the saturation phase.
As can be seen from Fig. 5a, the bar instability in the real Nbody 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 barlike perturbation which rotates with a welldefined 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 ω_{I} ≈ 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, B1, B1m, B4, B4m and B23. 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 B0 (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 B1m and B4m show a delay in bar formation, while in the fullyrescaled runs, B1, B4, and B23, the bar forms earlier compared to B0. The time difference in the bar formation Δ_{B} is given in Table 3.
Fig. 6 Bar formation in models with satellites: a) bar amplitudes A_{2}/A_{0} are given for B0, B1, B1m, B4, B4m and B23 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 massrescaled runs, the initial positions of the satellites are further away from the host galaxy than in the fullyrescaled runs, so the interaction with the disc peaks later in the mass rescaled runs. In particular, the minimum distance in the B1m run is at 0.96 Gyr, while in the B1 run it occurs at 0.85 Gyr. The delay Δ_{B} is larger in the B4m run where four satellites are used.
The fullyrescaled models show contrary effect on the bar formation, compared to the massrescaled simulations. In the B1 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 B4 run where the four largest satellites are used, the time difference is 300 Myr. A comparison with the B23 run shows that the most massive satellites determine mainly the time difference – it is only 10 Myr larger than in the B4 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: (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 inplane 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 flyby, 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 zeropressure discs, because a cold stellar disc is dynamically equivalent to a fluid disc with zero pressure. If the perturbation has mfold 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):
where B ≡ −κ^{2}/ 4Ω is the Oort constant; Σ_{0} is the disc surface density and Φ_{a} is the potential perturbation; (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 (21)with phase function F(R), where k = dF(R) / dR and  kR  ≫ 1. From Eq. (19)one has (22)that is, surface density of the wave is in phase with the radial velocity outside corotation, and in antiphase inside corotation: (23)The equation for the azimuthal velocity (Eq. (18)) simplifies to (24)and taking into account that Φ_{a} = −2πGΣ_{da}/  k , (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 B0 run (no satellites) at the time of the encounter in B1 (the pink diamond shows the position of the primary satellite in B1 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, (26)that is, (27)meaning that radial velocity lags behind the surface density by π/ 4 (m = 2). For the azimuthal velocity, (28)that is, (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_{2}(R,t) /A_{0}(R,t), (30)and velocity components, V_{k,2}(R,t) /A_{0}(R,t), (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 B0 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.
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 zoomin of the early growth phase for the B0, B1, and B1m runs. In particular, we stress the following features.
In the fullyrescaled B1 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 massrescaled model B1m, 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 B0 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 barmode 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 B0, as well as for one satellite runs B1 (fullyrescaled) and B1m (massrescaled). The lower left panel contains the map for notyetmentioned B1110 run, which differs from B1 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 massrescaled B1m run (see Fig. 3).
Fig. 11 Surface density maps log  A_{2}/A_{0}  (left panels, red = high, blue = low amplitude); Re A_{2}/A_{0} (right panels, red = positive, blue = negative amplitude at the xaxis) for B0, B1, B1m, and B1110 runs (top to bottom), reflecting the interference between the barlike 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 B0 run, it is the interference between a small barlike 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 B1, B1m, B1110. In Fig. 14 one can see for case B1 the onset of the disc perturbation after the encounter (0.87 and 0.89 Gyr), which clearly has a barlike 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(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 timescale 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_{2}/A_{0}, that is, a joint m=2amplitude along the xaxis, for the same runs B0, B1, B1m, B1110 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,(33)meaning that R = R(k,ω) is a function of k; ω ≡ mΩ_{s} is a parameter. For very open spiral^{2}, 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Ω)^{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: (34)where k is the wavenumber and c is a sound speed, which we use interchangeably with the radial velocity dispersion. Thus (35)and substitution to the WKB relation gives (36)In the last expression, Q denotes Toomre stability parameter for fluid discs, (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 B1 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 B1 run is definitely stronger than its counterpart in the B0 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 B1, it also works for enhancement.
In the last two rows, we see the opposite behaviour. In B1m, the red bar spot immediately after t = 1 Gyr is definitely weaker than in B0 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 B0 panels between 1.05 and 1.3 Gyr, and 1.15 and 1.4 Gyr. They also came into play and decrease amplitudes in B1m and B1110 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 B1500 run, in which the fully scaled satellite approaches the host galaxy 500 Myr earlier than in the B1 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 B0 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 B1500 curve continues almost in pace with the B1 curve.
Fig. 13 Bar amplitudes for B0 and B1500 runs. 
Fig. 14 Onset of the largescale barlike perturbation from the satellite in the B1 run (time from top to bottom, Σ_{d}, V_{R} and V_{θ}from left to right). 
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 AqD2 is believed to be a fair representative of a Milky Waylike host. Figure 3 in Moetazedian & Just (2016) is very helpful for comparing the distribution of satellites from the AqD2 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^{9}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^{8}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 AqD2, 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^{9}M_{⊙} (Bekki & Stanimirović 2009) and galactocentric distance of 60 kpc (McConnachie 2012) have a corresponding analogue in the AqD2 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^{2}. For a Sgr or SMClike 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 Nbody cosmological simulations and observational studies are dominated by lower mass objects. In the Aquarius satellites distribution there exists a fair number of such lowmass satellites with M_{tid}< 10^{9}M_{⊙} and pericentres within the inner 10 kpc of the disc. Observations of the Local Group have shown the existence of socalled ultrafaint 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 ultrafaint dwarfs to 17 (DrlicaWagner 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 yeartwo quick release (Y2Q1) predicts a total of ~100 ultrafaint 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 ultrafaint 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 Waylike galaxy model is studied mainly using Nbody 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 efolding 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 efolding 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 razorthin 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 barlike 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 D2 cosmological simulation. All satellites with tidal mass M_{tid} above 10^{8}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^{9}M_{⊙}) and one run with all 23 satellites. Each satellite was represented by 50 000 particles. In all cases, the pattern speed and the efolding 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 (B1m, B1110), or advancement (B1), 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 B1 and B1110 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 ≳ t_{2}. If the encounter occurs well before t_{1}, the tidal wave decays with no effect. This is proved by our B1500 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 B0. During the period of vulnerability, t_{1} ≲ t ≲ 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 longlasting 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 cofunded 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 041043 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 155212387, 160200649, and by the Basic Research Program OFN15 “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
 Abraham, R. G., Merrifield, M. R., Ellis, R. S., Tanvir, N. R., & Brinchmann, J. 1999, MNRAS, 308, 569 [NASA ADS] [CrossRef] [Google Scholar]
 Athanassoula, E. 2005, MNRAS, 358, 1477 [NASA ADS] [CrossRef] [Google Scholar]
 Athanassoula, E., & Sellwood, J. A. 1986, MNRAS, 221, 213 [NASA ADS] [CrossRef] [Google Scholar]
 Baugh, C. M., Cole, S., & Frenk, C. S. 1996, MNRAS, 282, 27 [NASA ADS] [CrossRef] [Google Scholar]
 Bédorf, J., Gaburov, E., & Portegies, Z. S. 2012a, in Advances in Computational Astrophysics: Methods, Tools, and Outcome, eds. R. CapuzzoDolcetta, M. Limongi, & A. Tornambè, ASP Conf. Ser., 453, 325 [Google Scholar]
 Bédorf, J., Gaburov, E., & Portegies, Z. S. 2012b, Astrophysics Source Code Library [record ascl: 1212.001] [Google Scholar]
 Bekki, K., & Stanimirović, S. 2009, MNRAS, 395, 342 [NASA ADS] [CrossRef] [Google Scholar]
 Bertin, G. 2014, Dynamics of Galaxies (Cambridge: Cambridge University Press) [Google Scholar]
 Binney, J., & Tremaine, S. 2008, Galactic Dynamics (Princeton: Princeton University Press) [Google Scholar]
 BlandHawthorn, J., & Gerhard, O. 2016, ARA&A, 54, 1 [Google Scholar]
 Borne, K. D., Bushouse, H., Colina, L., et al. 1999, Ap&SS, 266, 137 [NASA ADS] [CrossRef] [Google Scholar]
 DrlicaWagner, A., Bechtol, K., Rykoff, E. S., et al. 2015, ApJ, 813, 109 [NASA ADS] [CrossRef] [Google Scholar]
 Dubinski, J., Gauthier, J.R., Widrow, L., & Nickerson, S. 2008, ASP Conf. Ser., 396, 321 [NASA ADS] [Google Scholar]
 Dubinski, J., Berentzen, I., & Shlosman, I. 2009, ApJ, 697, 293 [NASA ADS] [CrossRef] [Google Scholar]
 Elmegreen, D. M., Elmegreen, B. G., & Bellin, A. D. 1990, ApJ, 364, 415 [NASA ADS] [CrossRef] [Google Scholar]
 Gillessen, S., Eisenhauer, F., Trippe, S., et al. 2009, ApJ, 692, 1075 [NASA ADS] [CrossRef] [Google Scholar]
 HernándezToledo, H. M., AvilaReesev, V., Conselice, C. J., & Puerari, I. 2005, AJ, 129, 682 [NASA ADS] [CrossRef] [Google Scholar]
 Hohl, F. 1971, ApJ, 168, 343 [NASA ADS] [CrossRef] [Google Scholar]
 Kalnajs, A. J. 1971, ApJ, 166, 275 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Kalnajs, A. J. 1977, ApJ, 212, 637 [NASA ADS] [CrossRef] [Google Scholar]
 Kannan, R., Macciò, A. V., Fontanot, F., et al. 2015, MNRAS, 452, 434 [NASA ADS] [CrossRef] [Google Scholar]
 Kazantzidis, S., Bullock, J. S., Zentner, A. R., Kravtsov, A. V., & Moustakas, L. A. 2008, ApJ, 688, 254 [NASA ADS] [CrossRef] [Google Scholar]
 Kormendy, J. 1982, ApJ, 257, 75 [NASA ADS] [CrossRef] [Google Scholar]
 Kuijken, K., & Dubinski, J. 1995, MNRAS, 277, 1341 [CrossRef] [Google Scholar]
 Law, D. R., Johnston, K. V., & Majewski, S. R. 2005, ApJ, 619, 807 [NASA ADS] [CrossRef] [Google Scholar]
 Lora, V., Grebel, E. K., SánchezSalcedo, F. J., & Just, A. 2013, ApJ, 777, 65 [NASA ADS] [CrossRef] [Google Scholar]
 Lotz, J. M., Jonsson, P., Cox, T. J., & Primack, J. R. 2010, MNRAS, 404, 575 [NASA ADS] [CrossRef] [Google Scholar]
 McConnachie, A. W. 2012, AJ, 144, 36 [NASA ADS] [CrossRef] [Google Scholar]
 Martin, N. F., Geha, M., & Ibata, R. A. 2016, MNRAS, 458, 59 [Google Scholar]
 MartinezValpuesta, I., Shlosman, I., & Heller, C. 2006, ApJ, 637, 214 [NASA ADS] [CrossRef] [Google Scholar]
 Miller, R. H., Prendergast, K. H., & Quirk, W. I. 1970, ApJ, 161, 903 [NASA ADS] [CrossRef] [Google Scholar]
 Moetazedian, R., & Just, A. 2016, MNRAS, 459, 2905 [NASA ADS] [CrossRef] [Google Scholar]
 Monari, G., Famaey, B., Siebert, A., et al. 2017, MNRAS, 465, 1443 [NASA ADS] [CrossRef] [Google Scholar]
 Moore, B., Ghigna, S., & Governato, F. 1999, ApJ, 524, 19 [Google Scholar]
 Noguchi, M. 1987, MNRAS, 228, 635 [NASA ADS] [Google Scholar]
 Navarro, J. F., Frank, C. S., & White, S. D. M. 1997, ApJ, 490, 493 [NASA ADS] [CrossRef] [Google Scholar]
 Ostriker, J. P., & Peebles, P. J. E. 1973, ApJ, 186, 467 [NASA ADS] [CrossRef] [Google Scholar]
 Polyachenko, E. V. 2005, MNRAS, 357, 559 [NASA ADS] [CrossRef] [Google Scholar]
 Polyachenko, E. V. 2016, Balt. Astron., 25, 288 [NASA ADS] [Google Scholar]
 Polyachenko, E. V., Berczik, P., & Just, A. 2016a, MNRAS, 462, 3727 [NASA ADS] [CrossRef] [Google Scholar]
 Polyachenko, E. V., Berczik, P., & Just, A. 2016b, Balt. Astron., 25, 411 [NASA ADS] [Google Scholar]
 Polyachenko, V. L., Polyachenko, E. V., & Strelnikov, A. V. 1997, Astron. Lett., 23, 525 [NASA ADS] [Google Scholar]
 Reid, M. J., Menten, K. M., Brunthaler, A., Zheng, X. W., & Dame, T. M. 2014, ApJ, 783, 130 [Google Scholar]
 Simon, J. D., & Geha, M. 2007, ApJ, 670, 313 [NASA ADS] [CrossRef] [Google Scholar]
 Simon, J. D., Li, T. S., & DrlicaWagner, A. 2017, ApJ, 838, 11 [NASA ADS] [CrossRef] [Google Scholar]
 Springel, V., Wang, J., Vogelsberger, M., et al. 2008, MNRAS, 391, 1685 [NASA ADS] [CrossRef] [Google Scholar]
 Toomre, A. 1969, ApJ, 158, 899 [NASA ADS] [CrossRef] [Google Scholar]
 Toomre, A. 1981, in Structure and evolution of normal galaxies, ed. S. M. Fall, & D. LyndenBell (Cambridge Univ. Press.), 111 [Google Scholar]
 Van den Bergh, S., Abraham, R. G., Ellis, R. S., et al. 1996, AJ, 112, 359 [NASA ADS] [CrossRef] [Google Scholar]
 Vogelsberger, M., Genel, S., Springel, V., et al. 2014, MNRAS, 444, 1518 [NASA ADS] [CrossRef] [Google Scholar]
 Walker, I. R., Mihos, J. C., & Hernquist, L. 1996, ApJ, 460, 121 [NASA ADS] [CrossRef] [Google Scholar]
 Wegg, C., Gerhard, O., & Portail, M. 2015, MNRAS, 450, 4050 [NASA ADS] [CrossRef] [Google Scholar]
 Widrow, L. M., Pym, B., & Dubinski, J. 2008, ApJ, 679, 1239 [NASA ADS] [CrossRef] [Google Scholar]
 White, S. D. M., & Rees, M. J. 1978, MNRAS, 183, 341 [NASA ADS] [CrossRef] [Google Scholar]
 Wu, Y.T., Pfenniger, D., & Taam, R. E. 2016, ApJ, 830, 111 [NASA ADS] [CrossRef] [Google Scholar]
 Zinchenko, I. A., Berczik, P., Grebel, E. K., Pilyugin, L. S., & Just, A. 2015, ApJ, 806, 267 [NASA ADS] [CrossRef] [Google Scholar]
All Tables
Parameters of the bar (bar radius, pattern speed, corotation radius, ellipticity) in B0 run at different moments in time.
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 M_{tid}> 10^{8}M_{⊙} (black) and M_{tid}> 10^{9}M_{⊙} (blue) before applying the rescaling. The red arrow marks the position of the AqD2 encounter used in simulations with one satellite and it is the 5th encounter from the right handside of the upper cumulative line. The satellite with the highest impact (~ 0.4) has a mass 1.1 × 10^{8}M_{⊙} 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 B0 run oriented along the xaxis 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 A_{2}/A_{0} and U_{Rθ,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 A_{2}/A_{0} are given for B0, B1, B1m, B4, B4m and B23 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 B0 run (no satellites) at the time of the encounter in B1 (the pink diamond shows the position of the primary satellite in B1 for visualisation). 

In the text 
Fig. 9 Same as in Fig. 8, but at the time of the encounter in B1m. 

In the text 
Fig. 10 Bar amplitudes  A_{2}/A_{0}  for the B0, B1, B1m runs (selected zoom of Fig. 6). 

In the text 
Fig. 11 Surface density maps log  A_{2}/A_{0}  (left panels, red = high, blue = low amplitude); Re A_{2}/A_{0} (right panels, red = positive, blue = negative amplitude at the xaxis) for B0, B1, B1m, and B1110 runs (top to bottom), reflecting the interference between the barlike 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 B0 and B1500 runs. 

In the text 
Fig. 14 Onset of the largescale barlike perturbation from the satellite in the B1 run (time from top to bottom, Σ_{d}, V_{R} and V_{θ}from left to right). 

In the text 
Fig. 15 As in Fig. 11 for B0 and B1500 runs for the period 0.2 < t < 1.4 Gyr. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.